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Abstract 


This  report  documents  codes  for  the  numerical  solution  of  control 
and  optimal  control  problems  for  diffusion  or  reflected  diffusion  mod¬ 
els  in  dimensions  two  to  four  and  for  continuous  time  Markov  chain 
control  problems  where  the  state  space  of  the  chain  is  a  grid  in  such 
a  Euclidean  space.  The  control  appears  linearly  in  the  dynamics  and 
cost  function  but  otherwise  the  process  and  cost  function  are  general. 
The  underlying  numerical  methods  use  efficient  forms  of  the  approx¬ 
imation  in  policy  space  and  multigrid  type  methods,  based  on  the 
Markov  chain  approximation  method  of  [7]. 
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1  Introduction 


Numerical  methods  of  stochastic  control  have  been  successfully  applied  in 
many  areas.  The  further  use  of  these, techniques  requires  the  availability  of 
easy  to  use  and  flexible  codes,  which  can  be  adapted  to  a  variety  of  basic  prob¬ 
lems.  This  report  documents  a  set  of  codes  which  were  developed  originally 
for  the  numerical  study  of  various  problems  in  modern  telecommunications. 
Those  applications  were  quite  successful  [5,  8].  These  references  discussed 
a  versatile  numerical  method  for  getting  good  approximations  to  queueing 
and  multiplexing  systems  and  getting  optimal  controls,  and  various  simple 
approximations  to  them,  as  well  as  for  experimenting  with  approximation 
schemes  when  the  systems  data  are  not  well  known.  Much  supporting  nu¬ 
merical  data  was  given  and  it  was  shown  that  the  technique  can  be  of  great 
help  in  the  design  and  analysis  of  systems. 

The  codes  developed  for  those  applications  are  quite  general  in  nature 
and  can  be  applied  to  a  great  variety  of  stochastic  control  problems  with 
diffusion  or  reflected  diffusion  type  models.  They  can  also  be  used  for  the 
solution  of  control  problems  for  continuous  time  Markov  chain  models  whose 
state  spaces  are  regular  grids  in  Euclidean  space. 

This  report  is  a  documentation  and  user’s  guide  to  the  codes  for  work¬ 
stations.  The  user  must,  of  course,  supply  the  data  for  the  system  and  cost 
function  of  interest.  The  format  will  be  discussed  in  detail  in  the  report 
and  will  be  illustrated  by  concrete  examples.  On  the  boundary  of  the  state 
space,  the  process  can  be  either  reflecting  or  absorbing,  as  desired.  The  pro¬ 
gram  can  be  used  in  many  ways,  as  described  in  the  sequel.  The  choices 
are  indicated  to  the  program  bj'^  an  options  parameter,  which  is  entered  at 
run  time.  The  options  include  various  choices  over  the  input  and  output. 
In  many  problems  the  cost  function  of  interest  is  the  sum  of  several  terms. 
While  one  might  want  the  optimum  policy  and  cost  for  the  total  cost,  the 
individual  components  of  the  cost  are  often  of  great  interest  as  well.  These 
can  be  evaluated  under  the  optimal  control.  Additionally,  one  has  the  choice 
of  not  optimizing,  but  merely  evaluating  the  cost  (and  its  components)  under 
a  given  control.  Threshold  controls  are  of  great  interest,  and  it  is  easy  to 
use  the  code  to  evaluate  performance  under  such  fixed  controls.  In  addition, 
the  numerical  solutions  and  controls  can  be  saved  for  postprocessing  analysis 
such  as  control  curve  plotting.  More  details  are  given  in  Section  3.6  which 
lists  the  run  time  options.  The  flexibility  offered  by  this  method  allows  a 
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wide  range  of  problems  to  be  computed  with  hopefully  minimal  effort.  High 
performance  machines  such  as  the  Cray  C90  require  (limited)  modifications 
to  achieve  effective  vector  performance  and  reduced  execution  times. 

Efficient  numerical  solution  techniques  are  used  to  get  the  various  values 
of  interest  and  the  optimal  controls.  The  original  process  and  cost  function 
are  approximated  by  the  versatile  Markov  chain  approximation  method.  The 
basic  numerical  scheme  then  uses  the  approximation  in  policy  space  method 
for  solving  these  approximating  control  problems.  This  approximation  in 
policy  space  method  generates  a  minimizing  sequence  of  policies.  For  each 
such  policy,  there  is  a  linear  system  of  equations  which  must  be  solved,  and 
this  is  the  core  of  the  computational  burden.  This  system  of  equations  rep¬ 
resents  the  cost  function  for  a  Markov  chain  control  problem.  Experience 
has  shown  that  overrelaxed  Gauss- Seidel  iterations  combined  with  multigrid 
methods  work  very  well,  and  (generally)  significantly  reduce  overall  com¬ 
putation  time.  Multigrid  use  without  overrelaxation  also  works  well.  The 
software  allows  the  user  to  tailor  the  computations  by  specifying  the  number 
of  multigrid  sublevels  and  the  overrelaxation  parameter  for  each  level. 

The  basic  diffusion  model  is  discussed  in  Section  2.  The  description  cen¬ 
ters  around  the  process  with  reflecting  boundaries  since  it  is  more  compli¬ 
cated  than  the  absorbing  boundary  case.  For  the  reflecting  boundary  case, 
the  cost  can  be  either  of  the  ergodic  or  discounted  type.  The  details  con¬ 
cerning  the  structure  of  the  state  space  and  the  basic  numerical  methods  are 
similar  to  those  for  the  absorbing  boundary  problem,  and  the  few  changes 
and  simplifications  (adding  the  absorbing  boundary  cost  and  dropping  the 
reflection  terms  and  associated  costs)  are  mentioned  at  the  end  of  the  sec¬ 
tion.  The  state  space  is  a  hyperrectangle  in  all  cases.  Section  3  is  concerned 
with  the  definitions  needed  by  the  program,  and  the  input  data  formats  and 
the  files  which  the  user  provides.  The  codes  are  very  flexible,  and  the  struc¬ 
ture  was  developed  so  that  many  types  of  problems  could  be  accommodated. 
Section  3  also  explains  how  to  compile  and  run  the  program,  and  how  to 
select  the  various  options,  as  well  as  the  possible  outputs.  A  detailed  exam¬ 
ple  is  provided  to  illustrate  the  procedure.  Sections  4  and  5  contain  other 
illustrative  examples.  The  above  examples  are  for  the  reflecting  boundary 
case.  An  example  for  the  absorbing  boundary  case  is  provided  in  Section  6. 
Section  7  deals  with  the  continuous  time  Markov  chain  control  problem.  The 
procedure  is  essentially  the  same  as  for  the  diffusion  model,  except  that  the 
controlled  transition  rates  must  be  defined.  Section  8  contain  an  example  of 
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a  continuous  Mai'kov  chain  problem  with  controlled  transition  rates. 


2  Process  Model  and  Cost  Function  Descrip¬ 
tions 

The  codes  can  be  used  for  the  numerical  solution  of  optimal  control  problems 
where  the  underlying  model  is  of  the  diffusion  or  reflected  diffusion  type  in 
dimensions  two  to  four  and  for  the  calculation  of  costs  or  cost  functions  associ¬ 
ated  with  given  controls.  They  can  also  be  used  for  the  Markov  chain  models 
where  the  state  space  is  a  regular  grid  in  a  Euclidean  space  of  dimensions 
two  to  four.  In  the  next  subsection,  we  describe  the  diffusion  models.  In  all 
cases,  the  state  space  is  a  hyperrectangle  G  in  Euclidean  d— space,  d  —  2,3, 4. 
The  particular  control  problems  which  originally  motivated  the  development 
of  the  codes  arose  in  queueing  or  telecommunications  situations,  where  the 
natural  state  space  was  often  a  hyperrectangle.  In  these  applications,  the 
states  frequently  correspond  to  (scaled)  buffer  or  queue  occupancies,  and 
hence  are  non-negative.  There  are  also  physical  upper  bounds,  which  yield 
the  upper  bounds  of  the  confining  rectangle.  The  reflection  directions  on  the 
boundarj^  were  determined  by  the  physics  of  the  flow  within  the  network  and 
were  constant  on  each  face  of  the  hyperrectangle.  The  present  code  keeps 
the  general  structure,  hence  the  state  is  confined  to  some  hyperrectangle. 

To  do  numerical  work,  one  needs  to  work  in  a  bounded  state  space.  If  the 
basic  state  variables  in  a  model  do  not  have  natural  finite  bounds,  then  they 
generally  need  to  be  bounded  in  some  fashion,  so  that  “finite”  procedures 
can  be  used.  The  most  common  approach  is  the  artificial  (upper  and  lower) 
truncation  of  each  state  at  some  appropriately  large  values.  To  prevent  the 
process  form  e.xceeding  these  values,  one  uses  either  a  reflection  or  absorption 
on  the  boundary.  One  tries  to  select  the  boundary  behavior  such  that  it  does 
not  seriously  affect  the  numerical  data  of  most  importance.  Thus,  even  for 
general  numerical  problems,  one  might  still  wish  to  work  with  a  state  space 
which  is  a  hyperrectangle,  with  appropriate  boundary  conditions  imposed. 

The  hyperrectangle  state  space  is  defined  by 

G={x:X-  <  Xi  <  X+}, 

where  X^,X~  are  real  numbers.  The  cost  function  for  the  reflecting  or 
absorbing  boundary  diffusion  model  is  described  below  in  the  next  subsection. 
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2.1  Diffusion  Models:  Reflecting  Boundaries 


Structure  of  the  process  model.  The  drift  term  in  the  diffusion  can  have 
an  arbitrary  dependence  on  the  state  x,  but  it  is  assumed  to  be  linear  in  the 
control,  and  the  covariance  matrix  is  a  constant..  For  the  reflecting  bound¬ 
ary  case,  the  process  is  reflected  “inwards”  when  it  tries  to  leave  G.  More 
particularly,  there  are  vectors  p,,  and  qi,  which  are  the  reflection  directions 
on  the  surfaces  xi  =  X~  and  Xi  =  A’/",  resp.  These  are  not  necessarily  unit 
vectors,  and  we  will  now  explain  how  they  are  specified.  Define  the  matrices 
P  =  {Pi'  •  •  •  iPd},  Q  —  {?!,  • .  ■ ,  9d},  where  the  pi,  qi  are  the  columns.  Refer 
to  Figure  1  which  illustrates  a  two  dimensional  problem.  Suppose  that  we 
are  sitting  at  a  point  x  on  a  face  of  G  and  move  out  of  G  one  “unit”  in  the 
direction  orthogonal  to  the  surface,  and  denote  the  new  point  by  y.  The 
Pi^qi  are  the  vectors  needed  to  return  us  to  the  face  (hyperplane)  on  which 
X  lies  when  we  are  at  y  and  move  in  the  desired  reflection  direction.  The 
reflection  direction  depends  on  the  face  on  which  x  lies,  but  not  otherwise  on 
X,  i.e.,  they  are  constant  on  each  face.  This  is  the  typical  setup  in  problems 
which  arise  in  queueing  theory.  For  one  example,  let  d  =  3,  and  suppose  that 
the  face  is  defined  by  x\  =  .  Then  the  direction  vector  qi  takes  the  form 

qi  =  (  —  1,912,913)5  where  y  +  qi  is  on  the  hyperplane  defined  by  Xi  =  X^ . 
Note  that  the  first  component  must  be  negative.  If  the  hyperplane  is  defined 
by  Xi  =  A'f ,  then  we  have  x  =  y  +pi,  where  pi  =  (l,pi2,Pi3)-  Note  that  the 
first  component  must  be  positive.  With  this  format,  we  must  always  have 
l9ol  <  IP*jl  <  1- 

The  general  form  of  the  diffusion  process  can  be  written  as 

dx  —  fe(x,  u{x))di  +  crdW  +  PdL  +  QdU.  (2.1) 

The  term  b{x,u)  is  referred  to  as  the  “drift  function.”  The  drift  is  linear 
in  u  :  b{x,u)  =  b{x)  +  Ku,  where  K  is  a  matrix  of  real  numbers,  b{-)  is 
continuous  and  u(x)  is  a  feedback  control  with  M  real  valued  components 
u{x)  =  (ui(x), . . . ,  t/jvf(x)),  where  0  <  M  <  4.  There  are  Ui  such  that 
0  <  Uj(x)  <  Ui-  The  W{-)  is  a  standard  Wiener  process  (covariance  matrix 
being  the  identity).  The  covariance  matrix  E  =  aa'  is  a  constant.  The  L,  U 
are  vectors  whose  components  L,(-),  f7i(-),  ^  =  l,...,d,  are  non  decreasing 
real  valued  processes:  T(0)  =  U{0)  —  0,  and  Li{-)  can  increase  only  at  those 
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t  where  Xi{t)  =  X~ ,  and  Ui{-)  can  increase  only  at  those  t  where  Xi{t)  =  Xf . 
It  is  always  assumed  that  b{x)  and  the  direction  vectors  pi,qi  satisfy  what  is 
required  for  the  control  problem  to  be  well  defined  [7]. 


The  cost  function.  For  a  continuous  function  k[-),  define  the  “cost  rate” 
in  the  form 


N  M 

kix,  u)  =  ^  Ciki{x{t))  +  ^  CiUi{x{t)),  M  <  4.  (2.2) 

ii=:l  Z  =  1 

The  cost  function  can  be  of  either  the  ergodic  or  the  discounted  forms.  In 
program  usage,  the  c^,  Ci  are  combined  into  the  coefficients  c(l),  c(2),  etc.  For 
real  numbers  U^vi^  define  the  ergodic  cost  function 


k{x{t),  u{x{t)))dt  +  {hdLi{t)  +  VidUi{t)) 


(2.3) 


and  the  discounted  cost  function 


W{x,u)  r 

Jo 


k{x{t),  u{x{t)))dt  +  ^2  {hdLi{t)  -f  VidUi(t)) 


(2.4) 


where  the  initial  condition  is  a;(0)  =  x. 
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Computing  normalized  “components”  of  optimal  cost.  In  many  opti¬ 
mization  problems  the  cost  criterion  is  the  sum  of  several  components.  While 
the  minimum  value  is  important,  the  values  of  the  components  themselves 
might  also  be  important.  For  example,  in  a  queueing  problem  the  total  cost 
might  be  the  sum  of  its  components:  one  which  measures  delay  and  another 
which  measures  losses.  It  is  worthwhile  to  know  what  the  mean  delay  and 
the  mean  loss  are,  under  the  control  which  minimizes  the  sum  of  the  mean 
values. 

Suppose  that  the  optimal  control  and  cost  have  been  calculated.  One 
of  the  options  available  in  the  program  is  the  computation  of  the  values 
associated  with  the  components  of  the  costs:  the  values  of  (2.3)  or  (2.4) 
with  the  optimal  control  used  but  with  the  individual  components  of  the 
cost  replacing  the  bracketed  quantity  in  (2.3)  or  (2.4).  For  specificity,  let  us 
restrict  out  attention  to  (2.3).  If  computation  of  the  cost  associated  with 
the  components  is  requested  (see  below  for  instructions  on  how  to  make  this 
request  of  the  program),  the  cost  due  to  all  the  individual  Z,-,  Ui  in  (2.3)  with 
nonzero  coefficients  {li,Vi,  resp.)  will  be  computed  (but  with  the  coefficient 
li^Vi  set  equal  to  unity).  Similarly,  for  the  values  of  individual  ki  and  Uj- 
components,  for  all  i  such  that  c,  0  or  c,  ^  0.  The  total  number  of  nonzero 
values  of  c,-,Ci  can  be  at  most  3  *  dim,  where  dim  =  dimensionality  of  the 
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model. 

The  program  has  been  designed  to  be  flexible  and  easy  to  use.  In  order 
to  allow  for  reasonable  generality  without  excessive  code  complexity,  the 
program  has  been  organized  so  that  the  user  supplies  b{x,u)  and  k{x,u). 
This  is  done  in  a  simple  Fortran  code,  as  will  be  seen  in  the  examples  below. 
The  program  contains  numerous  options  from  which  the  user  can  choose  to 
tailor  the  computation  to  a  particular  problem. 

Comments  on  the  numerical  method.  For  the  numerical  solution,  the 
diffusion  model  (2.1)  with  the  given  cost  function  is  approximated  via  the 
Markov  chain  approximation  method.  See  [6,  7]  for  discussion  of  the  gen¬ 
eral  method.  The  latter  reference  contains  many  applications,  and  details 
concerning  approximations  and  numerical  methods.  At  present,  the  Markov 
chain  approximation  method  seems  to  be  the  method  of  choice  for  getting 
the  optimal  costs  and  controls  as  well  as  solving  for  costs  under  given  con¬ 
trols  for  general  stochastic  control  problems  with  diffusion  or  jump  diffusion 
type  models.  Discussions  of  the  actual  approximation  method  for  various 
multiplexing- type  systems  is  in  [8,  5]. 

The  idea  is  to  approximate  (2.1)  with  cost  (2.3)  or  (2.4)  by  a  suitable 
controlled  finite  state  Markov  chain  on  a  state  space  which  is  a  “discretiza¬ 
tion”  of  G.  With  suitable  choices  of  the  chain  and  associated  cost  function, 
the  value  of  the  associated  optimal  cost  approximates  that  of  the  original 
problem  arbitrarily  well.  (Similarly,  if  the  control  is  fixed.)  The  chain  is 
parameterized  by  a  parameter  h  (analogous  to  a  finite  difference  interval) 
such  that  as  /i  ^  0,  the  “local”  properties  of  the  chain  resemble  more  and 
more  closelj^  those  of  the  diffusion.  The  state  space  Gh  of  the  approximating 
chain  is  essentially  the  regular  /i— grid  on  G.  (Actually,  it  is  slightly  larger, 
the  total  number  of  grid  points  being 

n?,,  [(x+  -  x-]/h)  +  3] , 

with  two  of  the  extra  points  being  associated  with  the  “numerical  reflecting 
boundary.”  The  value  of  h  is  an  input  quantity,  to  be  supplied  by  the  user, 
as  described  below. 

Solution  techniques.  The  approximating  Markov  chain  control  problem 
is  solved  by  the  approximation  in  policy  space  method  [3,  7].  This  method 
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works  by  getting  a  sequence  of  control  policies  with  decreasing  cost,  until 
the  convergence  criterion  is  satisfied.  To  insure  against  pathologies,  the  user 
specifies  the  maximum  number  of  policy  updates  allowed.  For  each  control 
policy,  one  needs  to  get  an  approximate  solution  to  the  linear  equation  for 
the  cost  (discounted  cost  problem)  or  “relative  cost”  (ergodic  cost  problem). 
This  is  the  main  computational  work. 

The  user  can  specify  one  of  several  options  for  obtaining  these  approxi¬ 
mate  solutions.  The  simplest  method  is  by  use  of  a  Gauss-Seidel  relaxation. 
The  code  allows  the  use  of  overrelaxation,  which  generally  (but  not  always) 
provides  faster  convergence.  The  user  specifies  the  number  of  relaxations  per 
policy  update  as  well  as  the  overrelaxation  parameter.  To  use  this  direct 
relaxation,  set  the  “number  of  multigrid  sublevels”  input  parameter  to  zero. 
Good  values  of  the  overrelaxation  parameter  are  problem  dependent,  and  the 
user  must  do  a  little  experimentation.  The  best  values  tend  to  be  slightly  less 
than  the  value  at  which  the  algorithm  becomes  unstable.  For  the  problems 
in  [5],  we  generally  used  values  in  the  range  [1.1,1.25].  The  range  of  good 
values  tends  to  be  similar  for  problems  of  similar  structure,  so  good  values  for 
one  member  of  a  sequence  of  runs  for  closely  related  problems  will  generally 
be  good  for  the  other  members.  Assigning  values  of  1.0  to  the  overrelaxation 
parameters  results  in  simple  Gauss-Seidel  iterations. 

There  is  also  the  option  of  a  form  of  the  multigrid  method,  which  is 
generally  faster  [1,  2,  7],  and  which  we  use  whenever  possible.  The  number 
of  multigrid  sublevels  (i.e.,  grid  coarsenings)  must  be  specified,  up  to  a  user 
determined  maximum  value  supplied  at  compile  time.  The  present  default 
value  is  three.  See  Section  3.3  for  more  information  about  changing  the 
maximum  number  of  multigrid  sublevels.  If  “number  of  multigrid  levels” 
equals  zero,  then  only  a  Gauss-Seidel  procedure  (overrelaxed,  if  desired)  will 
be  used,  with  no  multigrid  sublevels.  Within  each  level  of  multigrid,  the 
code  uses  Gauss-Seidel  or  overrelaxed  Gauss-Seidel  smoothings.  The  user 
must  specify  the  number  of  relaxations  and  the  overrelaxation  parameter 
for  each  multigrid  level,  as  well  as  the  maximum  number  of  multigrid  cycles 
allowed.  The  code  does  a  policy  update  at  the  end  of  each  full  multigrid  cycle. 
Thus,  the  maximum  number  of  policy  updates  is  defined  by  the  maximum 
number  of  multigrid  cycles  allowed.  See  [4]  for  a  general  introduction  to  the 
multigrid  method,  and  [7]  for  a  brief  summary.  The  method  was  introduced 
into  the  stochastic  control  area  by  Akian  and  Quadrat  [1,  2].  Note  that  the 
use  of  multigrid  increases  the  required  memory,  up  to  at  most  twice  what  is 
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required  without  it. 

Program  execution  will  stop  when  either  the  maximum  number  of  control 
policy  updates  (or,  equivalently,  full  multigrid  cycles)  has  been  reached  or  the 
(sup)  norm  of  the  difference  between  successive  control  updates  has  reached 
its  desired  value,  whichever  comes  first. 

State  space  description.  The  hyperrectangle  G  is  described  in  the  follow¬ 
ing  way.  We  work  in  terms  of  an  “origin”  of  G,  and  the  discretization  level 
h.  This  approach  simplifies  the  code.  The  user  enters  the  number  of  grid  in¬ 
tervals  {N~ ,  N't]  for  the  negative  and  positive  directions  from  the  problem 
origin  for  each  of  the  coordinate  directions  ?,  as  well  as  the  discretization 
level  h,  and  the  location  Xq  of  the  origin  of  G.  Thus  Xf  =  Xo  -b  hNf  ■  One 
can  always  select  the  origin  to  be,  say,  the  “lower  left  hand  corner”  of  G,  in 
which  case  N~  =  0,  alH. 

Some  care  must  be  exercised  in  choosing  Xo^Nf ,  when  there  are  multi¬ 
grid  sublevels.  Each  successive  sublevel  of  the  multigrid  procedure  uses  a 
state  space  or  grid  whose  spacing  is  twice  that  of  the  previous  higher  level. 
Thus,  if  m  denotes  the  number  of  multigrid  sublevels  desired,  the  user  should 
ascertain  that  the  +  N~]/[2''],  k  =  l,...,miov  each  dimension  i  are  all 
integers.  If  more  multigrid  sublevels  are  specified  than  are  computationally 
consistent,  then  the  program  will  reduce  the  number  of  levels  which  are  to 
be  used  such  that  consistency  is  maintained.  Without  multigrid  usage,  the 
state  space  may  be  described  using  an  odd  number  of  mesh  intervals  in  any 
or  all  the  dimensional  discretizations  for  evaluation  by  overrelaxation. 

Threshold  and  arbitrary  controls.  It  is  often  desired  to  compute  costs 
associated  with  a  given  control.  The  program  can  be  used  for  this,  since 
the  user  can  input  any  desired  control  via  an  appropriately  formatted  data 
file.  This  control  can  be  used  in  two  ways.  It  can  serve  as  the  initial  control 
of  the  policy  updating  routine.  If  a  good  initial  controls  were  available, 
say  from  a  previous  run  for  related  problem,  this  will  save  computational 
time.  Alternatively,  one  can  chose  to  have  the  entered  control  remain  fixed 
throughout  the  computation,  in  which  case  the  program  will  compute  the 
values  of  the  selected  costs  under  that  control.  Input  control  is  enabled 
using  option  1  as  described  in  subsection  3.6  and  a  brief  description  of  the 
required  format  for  the  control  input  file  is  given  in  subsection  3.7. 

In  many  applications,  a  threshold  type  of  control  is  of  interest  due  to  the 
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simplicity  of  its  implementation.  To  facilitate  the  computation  of  costs  under 
a  threshold  control,  one  need  not  write  a  file  containing  all  the  control  values. 
All  that  one  needs  to  do  is  to  select  the  threshold  control  option  and  put  in 
the  appropriate  threshold  parameters  with  the  other  data.  A  description 
of  the  threshold  control  and  its  parametrization  is  given  under  option  16  in 
subsection  3.6. 


The  values  of  the  auxiliary  or  “relative  cost  function  for  the  er- 
godic  cost  problem.  With  the  Markov  chain  approximation  method,  the 
actual  functions  which  are  computed  are  costs  for  an  approximating  Markov 
chain  problem  (chosen  by  the  program,  using  methods  from  [7]).  In  order 
to  describe  the  quantities  computed,  it  is  useful  to  look  briefly  and  formally 
at  the  Bellman  equations  for  the  control  of  a  Markov  chain  with  an  ergodic 
cost  function.  Let  k{x,u)  be  a  cost  realized  by  the  chain  when  in  state  x  and 
control  u{x)  is  used.  For  the  ergodic  cost  problem,  the  Bellman  equation  is 
[3,  7,  9] 


V{x)  =  min 

a 


'^p{x,y\a)V{y)  +  k{x,a) -'y  , 
-  y 


(2.5) 


where  a  varies  over  the  set  of  allowed  control  values,  and  x,y  vary  over  the 
points  in  the  state  space.  7  is  the  optimal  average  cost  per  unit  time.  The 
function  V{x)  is  the  “relative  cost.”  For  a  fixed  control  u(-),  the  equation 
for  the  cost  and  relative  cost  W{x,u)  is 


lT(a;,u) 


^p(a;,y|u(a;))lT(t/,u)  +  A;(a;,u(a:))  -7(u)  , 
-  y 


(2.6) 


where  7(u)  is  the  average  cost  per  unit  time.  Most  of  the  computational  work 
consists  in  solving  a  sequence  of  equations  similar  to  (2.6).  The  solutions 
F(-)  or  W{-,u)  are  not  unique,  since  if  any  solution  is  modified  by  adding  a 
constant,  then  that  new  value  will  also  solve  (2.5)  or  (2.6),  resp. 

For  the  optimization  problem,  the  program  computes  7,  the  numerical 
approximation  to  the  optimal  average  cost  for  the  original  problem  (2.1), 
(2.3).  One  selects  a  “centering  point,  called  Xc,  discussed  in  more  detail 
below,  and  the  program  computes  the  relative  cost  F(a:)  —  V{Xc).  A  value  of 
Xc  must  always  be  chosen  for  the  ergodic  cost  problem.  The  values  V{x)  — 
V{Xc)  can  be  saved  in  a  file  (option  1024,  where  the  output  file  is  called 
value.data)  as  can  the  values  of  the  computed  controls  (options  1  and  512), 
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as  discussed  in  the  program  options  Section  3.6.  Analogous  comments  hold 
for  the  computation  of  'y{xi),W{‘,u)  for  the  fixed  control  case. 

The  centering  point  Xc-  For  the  ergodic  cost  problem,  the  matrix  P{u)  = 
{p{x,y\u{x)),x,y]  used  in  (2.5)  and  (2.6)  is  not  a  contraction  and  so  the 
equation  is  solved  slightly  indirectly.  The  algorithm  which  is  used  takes  ad¬ 
vantage  of  the  non-uniqueness  of  the  solutions  V{x)  and  W{x^u)  in  (2.5)  and 
(2.6),  resp.,  to  transform  the  equations  into  “contractions”  by  the  selection 
of  a  “centering  point”  Xc-  See  [7,  Chapter  7]  for  a  detailed  description  of 
the  method.  The  selection  of  Xc  is  critical  for  good  numerical  behavior  and 
some  experimentation  might  be  required.  Good  locations  for  one  problem 
are  typically  good  for  similar  problems.  Poor  choices  lead  to  poor  numerical 
behavior  (slow  convergence),  and  possibly  even  non-convergence.  The  value 
should  be  a  point  that  is  (loosely  speaking)  “as  recurrent  as  possible”  un¬ 
der  what  is  expected  to  be  the  optimal  control.  For  example,  it  should  not 
be  in  a  corner  which  is  visited  infrequently  relative  to  other  points.  Often, 
the  problem  has  an  inherent  stability.  For  example,  see  (3.1),  where  a  good 
choice  for  Xc  =  {Xc.i^Xca)  is  Xc,2  =  0  and  Xc,i  being  either  zero  or  slightly 
positive.  If  there  is  a  control  or  other  dynamical  force  which  pushes  the  path 
strongly  out  of  a  region,  then  AT  should  not  be  located  in  that  region. 

The  discounted  cost  problem  (2.4)  with  reflecting  boundaries.  For 

the  discounted  cost  problem,  one  may  specify  a  centering  point  Xc  or  omit  its 
usage.  This  is  a  run  time  option.  When  the  centering  point  option  is  used, 
the  numerical  algorithm  utilized  is  similar  to  the  one  used  for  the  ergodic  cost 
problem  discussed  above.  Care  must  be  exercised  with  the  choice  of  Xc  since 
a  good  selection  for  it  will  result  in  faster  convergence.  In  addition  to  the  use 
of  Xc  in  the  algorithm,  the  program  will  compute  V{Xc)  and  V(x)  —  V(Ac), 
with  the  latter  being  saved  in  a  file  if  desired.  The  output  was  chosen  in  this 
way  since  one  often  uses  small  discount  factors  which  lead  to  large  values 
of  the  costs,  but  moderate  values  for  the  differences  between  the  costs  at 
different  points.  It  seemed  preferable  to  center  the  values  and  print  the  large 
value  for  a  single  point  (namely,  Xc)  only. 

When  the  centering  point  option  is  not  used  (i.e.,  the  centering  point  is 
not  to  be  used  in  the  numerical  algorithm),  the  user  must  specify  a  sample 
point  As  at  input  time  in  lieu  of  the  centering  point.  This  is  for  “output” 
purposes  only.  The  program  outputs  the  value  C(As)  and  the  user  can  also 
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save  V(ic)  —  V(Xs)  in  a  file  (option  1024),  if  desired.  The  motivation  for 
computing  the  differences  of  V{x)  —  V(A^s)  is  the  same  as  stated  above. 


2.2  The  Absorbing  Boundary  Problem 

For  the  absorbing  boundary  problem  and  a  given  continuous  function  A:(-), 
the  “cost  rate”  is  defined  in  the  form  of  (2.2).  The  cost  function  is  then 
defined  as 

W{x,u)  =  E  f  [A;(a;(s),  u(a;(s)))ds]  +  Ee~^'^g{xr)  (2.7) 

JO 

where  r  =  min{f :  x{t)  ^  G}  and  initial  condition  a;(0)  =  x. 


3  Creating  an  Executable  Program 

3.1  Program  availability 

The  source  code  is  available  on  the  World  Wide  Web  from  the  Lefschetz 
Center  for  Dynamical  Systems  of  Division  of  Applied  Mathematics  at  Brown 
University’s  home  page  (http://www.dam.brown.edu/lcds.html).  Viewers 
should  select  the  “Software”  link  to  access  the  software  and  documentation. 
The  source  code  routines  and  example  include  files  have  been  bundled  to¬ 
gether  with  the  Unix  “shar”  command  into  a  single  Unix  shell  archive  file  for 
distribution.  Clicking  on  the  highlighted  “stochastic  control  software”  link 
will  initiate  a  file  transfer  of  the  bundled  source  code  to  the  user. 

Once  the  file  scmodels.shrl  has  been  received  by  the  user,  the  individual 
source  files  must  be  extracted  from  it.  The  command 

sh  scmodels.shrl 

is  the  means  by  which  the  files  are  extracted  from  the  archive  file  scmod¬ 
els.shrl.  This  command  will  create  a  directory  called  “sc_code”  into  which 
the  individual  source  code  files  will  be  placed.  Several  subdirectories  will 
also  be  created  which  contain  the  necessary  example  include  files  needed  to 
construct  the  model  programs  described  in  this  document.  Any  or  all  of 
the  example  programs  may  be  compiled  using  the  identical  makefile  sup¬ 
plied  within  each  of  the  example  subdirectories.  The  makefile  provided  with 
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the  source  code  has  been  configured  for  Sun  workstations  and  may  require 
some  editing  for  use  on  other  machines.  Reading  the  sample  include  files 
described  in  this  report  is  probably  the  best  method  to  familiarize  oneself 
with  the  style  and  method  of  utilizing. the  software.  The  sample  include  files 
can  be  copied  and/or  edited  to  construct  an  executable  program  for  the  the 
user-specific  model.  Comments,  bugs,  and  suggestions  should  be  directed  to 
djarvis@dam.brown.edu. 

3.2  How  to  Define  the  Model 

We  next  describe  the  method  of  putting  the  actual  model  information  into 
the  general  program,  and  then  show  how  to  run  the  program  and  select  from 
among  the  many  options.  In  order  to  provide  maximum  flexibility  without 
an  excessively  complex  code,  the  user  must  supply  certain  information  in  a 
fixed  (but  simple  and  reasonable)  format.  The  specific  features  of  the  model, 
such  as  the  drift  b(x,  u)  and  the  k{x,  u)  part  of  the  cost  function,  are  supplied 
to  the  code  through  the  use  of  Fortran  “include”  files  which  minimize  coding 
work  while  allowing  a  reasonable  flexibility  of  the  general  computational 
model.  This  user-supplied  code  should  follow  traditional  Fortran  standards, 
especially  the  fixed  source  code  form  requirement  that  statements  must  be 
within  column  positions  7  through  72.  Although  free  format  source  code  is 
accepted  by  many  compilers,  this  model  code  uses  a  fixed  form  to  increase 
portability.  There  are  a  total  of  thirteen  “include”  files  for  this  software 
and  three  of  these,  user  .timing. h,  userJimevar.h,  and  user. timeout. h  are 
used  for  optional  program  performance  timing.  The  file  userjprobs.h  is  used 
for  the  control  problem  when  the  original  process  of  interest  is  a  controlled 
Markov  chain  and  not  (2.1).  See  Sections  7  and  8  for  an  explanation  of 
the  continuous  Markov  chain  problem  and  its  implementation.  Problems 
having  absorbing  boundary  conditions  require  the  use  of  user  .boundary  .h 
to  describe  the  boundary  function  g{x).  See  Section  6  for  details  and  an 
example.  Otherwise,  the  remaining  seven  “include”  files  (listed  first  below) 
are  used  to  completely  describe  a  model  of  type  (2.1)  in  software.  Since  most 
compilers  will  seek  to  incorporate  the  include  files  during  the  compilation 
process,  all  the  include  files  should  be  present  in  the  directory  even  if  the 
files  themselves  are  not  being  used  (i.e.,  are  empty)  for  the  model. 

Fortran  “include”  files 
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•  user-drift. h:  specifies  the  drift  function  b{x,u) 

•  user_cost.h:  specifies  cost  k{x,u) 

•  user_covar.h:  calculations  of  covariance  S  for  model,  if  any 

•  user-var.h:  user  introduced  variables  in  above  files,  if  any 

•  userJn.h:  input  statements  for  user-introduced  variables,  if  any 

•  user_out.h:  output  statements  for  user-introduced  variables,  if  desired 

•  userJnit.h:  calculated  constants  or  initializations  for  model,  if  any 

•  user-timing. h:  machine  dependent  Fortran  timing  function,  if  used 

•  user-timevar.h:  user  introduced  variable(s)  for  timing  routine,  if  used 

•  user-timeout. h:  output  statement(s)  for  timing  data,  if  used 

•  user-boundary.h;  absorbing  boundary  function  g{x),  if  used 

•  user_kbdout.h:  output  statements  for  keyboard  entered  data,  if  desired 

•  user-probs.h:  transition  probabilities  for  the  Markov  chain  model 

There  is  a  required  notation  for  certain  variables  used  in  the  “include”  files 
because  these  files  are  incorporated  directly  into  the  code  and  the  notation 
must  be  consistent.  The  list  below  describes  these  variables  with  a  general 
outline  of  their  basic  usage.  An  explicit  example  then  follows.  Note  that  all 
but  one  of  the  provided  variables  are  indexed  with  the  user  providing  proper 
values  for  the  maximum  values  of  the  indices  in  most  cases.  In  the  use 
of  the  control  and  cost  variables,  the  index  i  is  used  as  a  system 
index  and  it  must  appear,  and  it  must  not  be  modified  by  the 
user.  The  index  i  represents  the  vector  spatial  index  for  a  given  point  which 
the  code  varies  over  all  the  mesh  points.  The  software  uses  vector  indexing 
to  represent  the  dimensional  data  in  an  explicit  linear  fashion  in  order  to 
increase  portability  and  performance. 

Fortran  Variables  with  Fixed  Notation 
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•  x(j),  j  =  l,  ■ .  •,  dim:  the  components  of  the  state  vector  x. 

The  variables  x(j)  appear  only  on  the  right-hand  side  in  the  user  s 
Fortran  statement  in  explicit  form  as  x(l),  x(2),  etc. 

•  u(i,j),  j=l,  . . . ,  M:  M  =  number  of  controls,  M  <  4. 

The  j-th  control  is  denoted  by  u(i,j).  The  variable  u(i,j)  appears  only 
on  the  right-hand  side  of  statements  in  the  user^drift.h  file  with  the 
index  j  used.  Index  i  represents  the  vector  spatial  index  and  must 
appear  and  must  not  be  modified  by  the  user.  The  code  will  iterate 
over  this  general  spatial  index. 

•  k(Lm),  (1=1... ..dim  is  the  coordinate,  m=l,...,M  is  the  control  compo¬ 
nent):  control  coefficient  matrix. 

This  is  the  matrix  K  in  the  definition  of  b{x,u)  in  (2.1).  Any  input 
entry  can  be  either  positive  or  negative.  The  parameter  k(l,m)  appears 
only  on  the  right-hand  side  of  statements  in  the  user.drift.h  file.  But 
all  entries  k(l,m)  of  the  matrix  K  must  be  assigned  at  data  input  time, 
even  if  their  values  are  zero. 

•  drift. 

In  usage,  this  variable  will  be  embedded  within  a  dimensionally  depen¬ 
dent  if-statement  which  is  indexed  by  the  variable  j.  See  the  examples. 
This  variable  appears  only  on  the  left-hand  side  of  statements  in  the 
userjdrift.h  file. 

•  cov(l,m),  (l,m  =  l,...,dim):  symmetric  covariance  matrix  S. 
User-supplied  input  values  or  Fortran  statements  to  assign  covariance 
values  to  the  specified  model.  This  file  is  needed  only  if  the  user  wishes 
to  represent  the  covariance  in  terms  of  other  parameters.  This  file 
gives  the  formulas  relating  them.  If  the  covariance  input  values  are 
simply  a  set  of  numbers,  then  the  file  can  be  left  empty  and  the  values 
inserted  with  the  other  data  at  run  time.  See  the  example  below.  For 
uncorrelated  covariance  (S  is  diagonal)  only  the  diagonal  elements  need 
to  be  supplied.  In  the  case  where  the  covariance  is  not  diagonal,  only 
the  nonzero  upper  tridiagonal  values  are  to  be  supplied.  The  program 
will  directly  assign  the  lower  diagonal  covariance  matrix  elements  from 
the  specified  upper  triangular  elements  to  assure  that  the  matrix  is 
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symmetrical.  The  variables  cov(l,m)  appear  only  on  the  left-hand  side 
of  statements  in  the  user.covar.h  file  and  all  indices  must  be  specified. 

•  c(l),  1=1,  . . .,  3*dim:  cost  coefficients. 

This  is  the  set  of  Ci,Ci  in  (2.2).  The  variables  c(l)  appear  only  on 
the  right-hand  side  of  statements  in  usev-cost.h  file  with  all  indices 
specified.  By  convention,  we  assume  that  the  number  of  control  cost 
coefficients  equals  the  number  of  controls  and  that  the  control  cost 
coefficients  appear  first  in  use. 

•  cost(i):  defines  the  cost  components  for  k{x,u). 

Variable  to  which  model  dependent  cost  formula  is  assigned.  Again, 
index  i  represents  the  vector  spatial  index  and  must  appear  but  not 
be  modified  by  the  user.  The  code  will  iterate  over  this  general  spa¬ 
tial  index.  The  variable  cost(i)  appears  only  on  the  left-hand  side  of 
statements  in  the  user. cost. h  file. 

•  value(i):  defines  the  absorbing  boundary  cost  at  a  calculated  vector 
boundary  index.  This  variable  appears  only  on  the  left-hand  side  of 
statements  in  the  user  .boundary  .h  file. 

•  xmin(j),  j=l,  ...,  dim;  the  minimum  boundary  value  for  dimension 
j  ,  as  evaluated  by  the  software.  This  variable  appears  only  on  the 
right-hand  side  of  statements  in  the  user  .boundary  .h  file,  if  needed. 

•  xmax(j),  j=l,  . . .,  dim:  the  maximum  boundary  value  for  dimension 
j,  as  evaluated  by  the  software.  This  variable  appears  only  on  the 
right-hand  side  of  statements  in  the  user  .boundary  .h  file,  if  needed. 

For  the  continuous  time  Markov  chain  model  described  in  Sections  7  and  8, 
the  user  will  also  have  to  define  the  controlled  transition  rates  r(a:,  ?/|u(a:)).  To 
facilitate  this  feature,  additional  system-defined  variables  are  made  available 
and  must  be  assigned  appropriate  values.  Refer  to  Sections  7  and  8  for  more 
information. 

Example.  In  order  to  illustrate  the  use  of  the  include  files  and  the  conven¬ 
tions  for  utilizing  both  system-provided  and  user-defined  variables,  consider 


18 


the  following  two  dimensional  problem,  which  originally  arose  as  a  heavy 
traffic  limit  to  a  multiplexer  problem  [8,  5]. 

dxi{t)  =  —  a)  dt  —  U)^[s)dt  -|-  dWi(t)  +  dLi{t)  —  dU\{t)  _ 

dx2{t)  =  —(A  +  n)x2{t)  dt  +  dW2{t)  +  dL2  —  dU2- 

The  covariance  matrix  is  diagonal  with  elements  cr(l,l)  =  0  and  cr(2,2)  = 
2  A  ///(A  +  p).  For  this  model,  the  drift  and  covariance  are  written  in  terms 
of  other  parameters:  A,a.  For  computational  purposes,  a  specific  value 

of  cr(2,2)  can  either  be  entered  directly  as  an  input  quantity  or  evaluated 
by  the  code  from  the  other  model  parameters  if  desired.  In  this  example, 
we  will  evaluate  the  covariance  values  via  the  “include”  file  to  illustrate  the 
appropriate  coding  procedure. 


Define  the  ergodic  cost  for  the  problem  as: 


7(u.)  =  lim  —E 
i\  /  T  T 


hUAT)  +  r  c(l)ut{s)ds  +  c(2)xt(s)di 


(3.2) 


The  “include”  and  data  files.  The  user  can  compile  one  of  two  versions  of 
the  program  for  entering  data.  In  the  “silent”  version,  no  message  prompts 
for  data  are  produced  by  the  program.  There,  the  user  prepares  an  input  file 
which  contains  all  of  the  required  data  in  the  order  specified  below.  In  the 
other  version,  called  the  “prompted”  style,  the  program  “asks”  the  user  to 
enter  the  necessary  system  variable  data  in  the  prescribed  order  by  producing 
data  message  prompts.  The  STDERR  descriptor  for  the  WRITE  statements 
in  the  file  userJn.h  below  is  used  for  the  “prompted”  input  version.  It  assures 
that  the  particular  message  prompts  supplied  by  the  user  will  be  written  to 
the  monitor  screen  for  entering  the  data  for  the  variables  of  the  problem  as 
the  program  sequences  through  its  input  statements. 

An  additional  feature  of  the  “prompted”  program  version  is  its  creation 
of  a  file  named  input. prompted  with  all  the  entered  system  data  values. 
A  complete  file  record  of  inputs  will  be  created  if  the  file  userJcbdout.h 
includes  the  necessary  WRITE  statements  using  the  KBD  output  descriptor 
to  write  the  input  data  for  user-introduced  variables  to  the  file.  The  file 
input. prompted  could  then  used  directly  as  an  input  file  for  a  subsequent  run 
of  the  program.  Each  execution  run  of  the  “prompted”  version  of  the  code 
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overwrites  an  existing  input. prompted  file  so  it  is  the  user’s  responsibility  to 
save  any  previous  input. prompted  files  of  interest. 

For  this  example,  we  will  assume  that  the  “prompted”  version  of  the  code 
will  be  compiled  by  the  user.  The  “include”  files  supplied  by  the  user  would 
be  as  written  follows: 

•  user_drift.h: 

if(  j  .eq.  1  )  then 

drift  =  nu='x(2)  -  a  -f  k(l,l)^u(i,l) 
else 

drift  =  —  (lambda  +  mu)  *  x(2) 
endif 

Comment:  The  j-th  component  of  the  spatial  variable  is  always  referred 
to  as  x(j).  The  variable  j  is  used  by  the  code  to  select  the  appropriate 
dimension  for  the  drift  evaluation  and  must  appear.  In  this  example, 
there  is  onlj^  one  control  and  and  it  is  denoted  by  u(i,l),  as  required. 
As  noted  above,  the  index  i  must  appear  as  given  above  for  proper 
indexing  of  the  control. 

•  user_cost.h: 

cost(i)  =  c(l)*u(i,l)  +  c(2)*x(l) 

Comment.  Note  that  the  parts  of  the  cost  dealing  with  the  reflection 
terms  L  and  U  are  not  included  here.  Since  their  structure  is  fixed, 
their  w^eights  are  entei-ed  with  the  other  parameters.  Again  note 
the  required  use  of  index  i  in  the  control  variable  u(i,l).  Otherwise, 
all  indices  appearing  in  the  cost  evaluation  have  been  specified  by  the 
user. 

•  user.covar.h: 

cov(2,2)  =  2.*lambda*mu  /  (lambda  +  mu) 


Comment:  In  this  example,  the  only  nonzero  covariance  value  is  given 
in  terms  of  other  parameters,  and  the  code  above  specifies  the  func¬ 
tional  dependence.  If  the  covariances  were  simply  a  set  of  given  real 
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numbers,  then  this  file  would  be  left  empty,  and  the  covariance  values 
would  be  entered  at  run  time  together  with  the  other  problem  data. 
Unless  the  user  uses  this  include  file  to  explicitly  evaluate  a  covariance 
value,  the  input  covariance  values  entered  at  run-time  are  used.  If  a 
particular  covariance  value  is  computed  in  this  include  file  and  is  also 
entered  as  data  at  run  time,  the  form  used  in  this  include  file  will  be 
the  one  used  by  the  program. 

•  user_var.h: 

real  lambda,  mu,  nu,  a 
common  lambda,  mu,  nu,  a 

Comment:  In  this  file  we  list  the  “user-supplied”  parameters  which 
were  introduced  in  the  above  files.  Values  for  these  variables  can  be 
entered  with  the  other  data  at  run  time  by  means  of  READ  statements 
or  can  be  assigned  by  the  user  in  the  user.init.h  file.  If  there  are  no  such 
parameters,  then  this  file  will  be  empty.  The  program  does  not  use  au¬ 
tomatic  typing  of  variables  so  all  introduced  variables  must  be  defined 
without  exception.  Compiler  complaints  about  undefined  variables  will 
follow  if  this  rule  is  ignored.  The  code  is  designed  to  use  64-bit  floating 
point  precision  in  its  calculations  for  improved  numerical  accuracy.  All 
real  (i.e.,  floating  point)  variables  for  a  given  model  should  be  typecast 
as  simple  Fortran  REAL  unless  some  specific  need  arises  for  another 
precision.  For  example,  the  variables  used  in  Sun  workstation  timing 
functions  require  .32-bit  precision.  See  the  comments  which  accompany 
the  userJimevar.h  file  usage  for  additional  example  information. 

By  using  the  “precision  neutral”  REAL  type  specification,  the  provided 
makefile  will  supply  the  appropriate  compiler  option  to  specify  that  all 
REAL  variables  are  to  be  64-bits  in  precision.  The  distinct  advantage 
of  linking  the  default  Fortran  REAL  variable  typing  with  the  64-bit 
precision  compiler  option  is  that  code  portability  is  enhanced  since  no 
assumption  is  made  with  respect  to  the  default  number  of  bits  used 
by  REAL  Fortran  variables  on  a  computational  platform.  This  re¬ 
solves  the  portability  and  maintenance  problems  that  using  DOUBLE 
PRECISION  declarations  pose  since  this  results  in  64  bits  on  many 
workstations  but  128  bits  on  high-performance  computers  such  as  the 
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Cray  C90. 

•  userJn.h: 

write(STDERR,*)  ’lambda  mu  nu  a  ?’ 
read(5,*)  lambda,  mu,  nu,  a 

Comment:  This  file  provides  the  necessary  statements  for  reading  data 
into  the  variables  which  the  user  has  introduced  for  the  model.  The 
inclusion  of  the  write  statement  is  for  the  case  where  the  user  desires 
prompted  keyboard  messages  for  data  input.  See  the  next  section  for 
more  on  the  use  of  STDERR  with  write  statements.  If  no  such  prompt¬ 
ing  is  desired,  omit  the  write  statement. 

•  user -Out.  h: 

write(6,'*')  ’lambda  =  ’,  lambda 
write(6,*)  ’  mu  =  ’,  mu 
write(6,*)  ’  nu  =  ’,  nu 
write(6,*)  ’  a  =  ’,  a 

Comment:  To  document  the  executed  problem,  this  file  allows  the  user 
to  provide  output  statements  to  print  variable  values  of  interest  which 
have  been  supplied  to  describe  the  model.  The  manner  and  style  in 
which  this  data  is  printed  is  left  to  the  user. 

•  userJnit.h: 

Comment.  This  illustrative  example  does  not  require  any  variable  ini¬ 
tializations  so  the  userJnit.h  file  would  be  left  empty.  To  demonstrate 
how  this  file  might  be  used,  suppose  that  we  have  a  model  with  a  cost 
defined  by  k{x,u)  =  (b  /  N)  g{x,u),  where  b,  Q,  and  N  are  in¬ 
put  parameters  that  might  vary  from  run  to  run  and  g{x,u)  is  some 
function.  We  may  reduce  the  execution  time  by  writing  k{x,u)  =  M 
g{x,u),  and  calculating  M  only  once: 

M  =  (b  /  N)  *  sqrt(Q) 
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where  M,  b,  Q,  and  N  have  been  declared  as  REAL  Fortran  variables 
in  the  file  user.vars.h.  Recall  that  REAL  variables  will  be  cast  with 
64-bit  representation  at  compile  time.  Additionally,  users  may  directly 
assign  values  to  their  model  dependent  variables  in  this  file. 

•  user_timing.h: 
time  =  dtime(t) 

Comment.  The  “dtime”  function  is  used  on  Sun  workstations  for  tim¬ 
ing  purposes.  The  variables  associated  with  the  function  as  well  as  the 
function  itself  must  defined  for  the  code.  For  computational  platforms 
other  than  Sun,  check  local  documentation  for  the  vendor-specific  tim¬ 
ing  function  call.  See  below  for  the  proper  declarations  of  the  timing 
variables  and  the  function  itself. 

•  user.timevar.h: 

reaR4  t(2),  time,  dtime 
common  /time/  t,  time 

Comment.  The  particular  timing  function  and  its  variables  are  defined 
for  the  code.  The  Sun  function  “dtime  ”  and  its  associated  variables 
require  32-bit  single  precision  real  variables.  To  enforce  this  precision 
we  use  the  REAL*4  declaration.  The  REAL*4  declaration  informs 
the  compiler  that  the  defined  variables  and  the  function  “dtime”  are  4 
bytes  (i.e...  32  bits)  in  precision.  Note  that  the  function  name  “dtime” 
is  omitted  from  the  Fortran  COMMON  statement  since  it  is  not  a 
variable. 

•  user.timeout.h: 

write(6,*)  ’Program  time  =  ’,  time(l)-(-time(2),  ’  secs’ 
write(6,*)  ’  user  time  =  ’,  time(l),  ’  secs’ 
write(6,*)  ’  system  time  =  ’,  time(2),  ’  secs’ 

Comment.  Print  the  results  of  the  timing  calls  for  the  program.  In  this 
example,  the  full  output  of  the  timing  calls  will  be  printed  as  specified 
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by  the  Sun  function  “dtime”.  Otherwise,  the  style  and  format  of  the 
output  is  left  to  the  individual  user. 

•  userJkbdout.h: 

write(KBD,*)  lambda,  mu,  nu,  a 


Comment:  This  file  provides  the  necessary  statements  for  writing  user- 
supplied  data  to  the  output  file  which  is  always  called  input. prompted 
and  which  is  created  when  the  user  compiles  the  code  for  prompted 
input.  The  data  should  be  written  in  the  same  order  as  it  is  read 
into  the  program.  Perhaps  the  simplest  way  to  generate  this  include 
file  is  by  changing  all  the  READ  statements  in  the  userJn.h  file  to 
WRITE  statements  which  use  the  KBD  descriptor.  Otherwise,  the 
style  and  format  of  the  output  is  left  to  the  individual  user.  Note 
that  the  descriptor  KBD  is  used  to  direct  the  output  to  the  default 
file  and  must  be  used  for  all  WRITE  statements.  If  the  program  has 
been  compiled  without  keyboard  prompting  then  omit  the  WRITE 
statements,  leaving  this  include  file  empty. 

Since  the  parameters  are  input  at  run  time,  the  proper  sequence  for  en¬ 
tering  the  data  will  be  deferred  until  we  have  discussed  the  compilation  pro¬ 
cedure.. 

3.3  Compiling  and  running  the  program 

The  UNIX  software  tool  “make”  is  used  to  compile  and  generate  the  ac¬ 
tual  executable  program.  The  “make”  command  is  directed  by  the  included 
makefile  in  compiling  and  linking  the  general  source  code  routines  into  the 
specified  program  with  as  little  user  intervention  as  possible.  Thus  the  user 
is  relieved  of  both  typing  and  remembering  the  compilation  commands  and 
options.  The  whole  compilation  and  link  process  will  be  apparent  to  the 
user  since  “make”  explicitly  describes  its  progress  as  it  executes.  Table  1  be¬ 
low  lists  the  appropriate  commands  to  issue  in  order  to  construct  the  listed 
program. 

The  “Input  mode”  column  in  table  1  represents  the  user’s  preferred  style 
of  entering  the  problem  data.  A  program  compiled  with  the  “Prompted” 
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Command 

Program 

Input  mode 

Dimension 

make  2d 

sctrl2d 

Silent 

Two 

make  3d 

sctrl3d 

Silent 

Three 

make  4d 

sctrl4d 

Silent 

Four 

make  2dask 

sctrl2dask 

Prompted 

Two 

make  3dask 

sctrl3dask 

Prompted 

Three 

make  4dask 

sctrl4dask 

Prompted 

Four 

Table  1:  Compilation  choices  for  the  control  code 

input  mode  will  issue  a  sequence  of  prompts,  at  each  of  which  the  user  will 
enter  the  appropriate  data  for  defined  system  inputs.  If  users  desire  program 
prompts  for  their  introduced  variables  then  users  must  provide  the  necessary 
WRITE  statments,  as  described  in  the  previous  section.  The  “Prompted” 
input  mode  will  also  save  all  system  and  introduced  variables  input  data  in 
a  file  named  input. prompted  for  future  reference.  (This  assumes  that  the 
user  has  supplied  the  necessary  WRITE  statements  in  the  file  userJthdout.h 
to  output  the  introduced  variable  data.)  In  the  “Silent”  input  mode,  the 
program  quietly  awaits  input  from  either  the  keyboard  or  a  file  and  no  de¬ 
fault  input  file  is  created.  For  either  program  version,  output  is  directed  (by 
default)  to  the  terminal  but  may  be  saved  by  redirecting  it  to  a  file  in  Unix 
fashion,  as  illustrated  below. 

Note:  Due  to  global  variable  dependencies  and  memory  layout,  it  is  neces¬ 
sary  to  utilize  the  “make  clean”  command  if  successive  runs  are  for  problems 
of  difi'erent  dimensions. 

Example.  As  an  example,  consider  the  case  where  the  user  desires  to  ex¬ 
periment  with  a  three  dimensional  model.  First,  let  the  input  mode  be 
“prompted.”  The  user  would  type  the  command 

make  3dask 

to  generate  the  prompted  three-dimensional  program.  [The  programs  for 
two-  and  four-dimensions  are  created  similarly.]  Once  the  “make  3dask” 
command  has  been  given,  the  source  code  will  be  compiled.  The  successful 
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compilation  and  linking  creates  the  program  named  “sctrl3dask,”  which  is 
the  one  to  be  executed  by  the  user.  Simply  running 

sctrlSdask 

will  cause  the  program  to  query  the  user  for  necessary  parameter  data,  save 
those  inputs  in  the  default  output  file  input. prompted  and  send  the  results  to 
the  terminal  screen.  [This  file  (perhaps  modified  by  the  user)  can  be  used  on 
a  subsequent  run  as  an  input  file,  to  save  time.]  Alternatively,  by  running 

sctrlSdask  >  user. output 

the  user  can  save  the  program  output  in  the  file  user. output  while  still  in¬ 
putting  data  at  the  program  prompts.  This  command  creates  (or  overwrites) 
the  output  file  user,  output.  The  output  file  contains  the  problem  parameters, 
stopping  data  (number  of  cycles,  etc.)  and  the  numerical  solutions.  The  user 
can  also  input  data  from  a  file  and  save  the  output  by  running 

sctrlSdask  <  input. data  >  user2. output 

thus  directing  data  from  the  file  input. data  to  the  “sctrl3dask”  program  and 
saving  (i.e.,  redirecting)  the  output  from  the  terminal  to  the  user-named  file 
user2. output.  This  is  the  preferred  method  of  entering  data,  especially  when 
using  the  “silent”  (i.e.,  unprompted)  program  versions,  since  it  is  the  quickest 
and  allows  the  user  to  verify  input  data  before  it  is  entered  and  used  by  a 
program.  Issuing  the  command 


make  clean 

command  removes  the  previously  created  “include”  files  (those  suffixed  with 
a  .h)  and  compiled  object  files  (suffixed  with  an  .o)  thereby  cleaning  out 
these  old  files  to  avoid  any  possible  global  memory  confusion,  and  saving 
disk  space  as  well. 

The  software  package  has  been  extensively  tested  within  its  Sun  SparclO 
workstation  development  environment  but  it  is  possible  that  compilation 
problems  will  appear  when  the  code  is  ported  to  other  machines.  Such  prob¬ 
lems  are  usually  the  result  of  local  system  features  such  as  the  Fortran  com¬ 
piler  name  and  library  function  names.  System  error  messages  can  be  very 
helpful  in  determining  the  nature  of  the  difficulty.  If  the  user  is  unable  to 
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resolve  a  compilation  problem,  it  is  suggested  that  help  be  secured  from  local 
support  personnel  or  more  experienced  Fortran  users. 

Important  default  settings.  There  are  three  important  default  settings 
at  compilation  time  which  the  user  should  carefully  review.  The  first  setting 
to  note  is  the  maximum  number  of  multigrid  sublevels  allowed.  The  present 
default  value  allows  three  multigrid  sublevels  (three  levels  of  multigrid  refine¬ 
ment).  Experience  has  demonstrated  that  two  or  three  multigrid  levels  tends 
to  result  in  the  shortest  execution  time  for  a  program.  This  default  value 
may  be  altered  by  changing  the  defined  GRIDS  value  in  the  included  makefile 
to  the  new  maximum  number  of  levels  desired.  Any  number  is  possible,  but 
the  number  is  an  upper  bound  to  what  can  be  entered  as  the  “number  of 
multigrid  sublevels”  (the  actual  numer  to  be  used  in  the  compiled  program) 
in  the  data  input. 

To  allow  prompted  input  messages  to  appear  for  keyboard  data  entry 
when  the  program  output  is  being  redirected  to  a  file,  the  program  writes 
these  messages  to  Fortran  standard  error  output  unit  (the  monitor).  On  the 
Sun  Sparc  10,  this  is  Fortran  unit  0.  Since  this  definition  of  the  standard 
error  output  unit  may  be  vendor  dependent,  this  value  may  be  altered  by 
changing  the  defined  ERR.OUTPUT  parameter  in  the  makefile.  The  value 
assigned  to  the  ERR_OUTPUT  parameter  will  be  inherited  by  the  STDERR 
output  descriptor  used  by  the  WRITE  statements  for  issuing  of  error  and 
warning  messages  as  well  as  interactive  keyboard  prompts,  if  used.  See  the 
vendor’s  Fortran  documentation  for  the  appropriate  local  standard  error  unit 
value.  As  a  practical  matter,  if  the  prompted  data  input  program  is  used 
and  the  messages  don’t  appear  on  the  monitor  screen,  then  the  value  for 
ERR-OUTPUT  is  probably  incorrect. 

The  third  default  setting  is  the  parameter  MAXJBITS.  This  value  rep¬ 
resents  the  highest  order  bit  value  for  machine  integer  representation  and  it 
is  hardware  dependent.  By  default,  the  parameter  MAXJBITS  is  set  to  31 
since  the  Sun  workstation  on  which  the  software  was  developed  and  tested 
uses  32  bit  integers.  For  most  workstations  this  value  of  MAXJBITS  will  be 
sufficient  since  the  32  bit  integer  representation  is  the  one  most  commonly 
used.  On  a  high  performance  machine  such  as  the  Cray  C90,  the  value  of 
MAXJBITS  would  be  set  to  63  since  the  standard  machine  word  length  of 
a  full  integer  is  64  bits. 
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3.4  Program  inputs 

As  noted  above,  the  input  data  can  be  entered  either  from  the  keyboard 
or  from  an  existing  file.  The  data  is  read  by  the  program  without  format 
restrictions  but  it  must  be  entered  in  the  order  specified  below,  and  using 
the  data  types  as  documented  below.  Failure  to  use  proper  data  types  may 
result  in  erroneous  results. 

The  stopping  criterion  tolerance  represents  the  maximum  absolute  dif¬ 
ference  for  the  finest  mesh  between  successive  policy  updates  or  between 
successive  full  multigrid  cycles  if  the  control  is  not  updated.  A  fairly  small 
value  for  this  parameter  should  be  chosen  (e.g.,  10“®  —  10“®). 

Values  for  the  “run  time  options”  input  parameter  are  described  further 
in  the  next  subsection.  The  following  input  description  is  for  the  reflecting 
boundary  case,  for  either  the  ergodic  or  the  discounted  cost  function.  Recall 
that  if  the  cost  is  of  the  discounted  type,  then  there  is  the  option  of  not  using 
the  centering  point  AT  in  the  numerical  algorithm.  Even  if  the  option  of  not 
using  the  centering  point  in  the  computation  is  elected,  then  an  input  value  of 
Xc  is  still  needed,  since  the  output  is  given  in  the  form  V{Xc).,  V(3:)  —  V(AT). 

General  input  data  order  and  data  types. 

Reflecting  Boundaries 

•  dimension  of  state  variable  (integer  =d) 

•  Number  of  mesh  intervals  [N^ ,  ),...,  )).  d  lines. (integers) 

•  Number  of  multigrid  sublevels  (integer) 

•  Number  of  controls  (M  <  4)  (integer) 

•  Non-diagonal  covariance  matrix?  (integer:  if  0  then  NO,  else  YES) 

•  h,  mesh  interval  width  (real) 

•  Origin  of  the  state  space  (real)  AT 

•  Value  for  run-time  options  (integer)  [See  Section  3.6  for  full  details] 

•  Centering  point  value  Xc  (real) 

•  Covariance  values:  upper  triangular  (real) 

d  lines:  ((7(1,1),...,  a{l,d)),  {a{2, 2), . . . ,  C7(2,  d)),. . . ,  {a{d,  d)) 

Note:  This  data  line  is  still  needed,  even  if  the  include  file  user^covar.h 
is  used.  Any  calculations  defined  in  user-covar.h  will  overwrite  the  cor¬ 
responding  values  entered  at  input  time. 

•  Drift  control  coefficient  matrix  K  (real:  d  x  (No.  of  controls)  matrix) 

This  entry  has  d  lines  containing  . . . ,  •  •  • ,  {kd,\,  ■  ■  ■ ,  kd,M)- 

•  Maximum  values  of  the  controls  ui, . . . ,  um  (real:  no.  of  controls) 
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•  Number  of  nonzero  cost  coefficients  (number  of  nonzero  c,,  c,  in  (2.2), 
called  c{i)  below)  (integer) 

•  Cost  coefficients  values  c(l), . . . ,  (real) 

•  User-supplied  model  inputs,  if  any  [In  the  example,  these  are  A,  /x,  n,  a] 
The  number  of  lines  used  is  specified  by  the  input  statement 

•  Cost  coefficients  for  underflow  /i, . . . ,  (d  real) 

•  Cost  coefficients  for  overflow  Ui, . , . ,  (d  real) 

•  Underflow  boundary  reflection  values,  (real:  d  times  d)  d  lines,  con¬ 
taining  pi ,....  pd 

•  Overflow  boundary  reflection  values,  (real:  d  times  d)  d  lines,  contain¬ 
ing 

•  Discount  factor  d  (real;  0.  =  no  discounting) 

•  Maximum  number  of  policy  updating  steps  (integer) 

•  Stopping  criterion  tolerance  (real) 

•  Number  of  relaxations  to  be  done  for  each  multigrid  level  (integers) 

•  Overrelaxation  factor  for  each  multigrid  level  (reals) 

•  Threshold  values  (real,  integer:  no.  of  controls)  These  inputs  are  read 
only  if  the  threshold  option  (option  16)  is  selected.  See  Section  3.6 
for  specifying  this  option.  There  are  M  lines,  corresponding  to  the  M 
controls.  With  the  threshold  control  option,  the  th  control  takes  the 
value  Ui  if  some  specified  state  component  equals  or  exceeds  a  given 
real  number.  Otherwise  the  i— th  control  takes  the  value  zero.  The 
z— th  line  is  (real,  integer),  the  integer  being  the  state  component,  and 
the  real  number  is  the  threshold  level.  If  some  control  is  never  to  be 
active,  set  the  threshold  level  large. 

3.5  Parameter  Inputs:  Example  1  Returned 

We  now  show  how  the  parameters  are  put  in,  either  via  a  file  or  via  the  mon¬ 
itor  and  keyboard.  For  any  prompted  program  version,  all  data  entered  from 
the  keyboard  will  be  written  to  a  default  output  file  named  input. "prompted 
which  is  created  by  the  program. 

The  first  six  inputs,  which  are  dimension,  number  of  mesh  intervals,  num¬ 
ber  of  multigrid  sublevels,  number  of  controls,  and  the  indicator  that  there  is 
a  non-diagonal  covariance  matrix  allow  the  program  to  determine  the  amount 
of  memory  which  is  to  be  allocated  for  the  problem.  This  tailoring  of  memory 
usage  for  a  problem  allows  efficient  memory  management  for  larger  memory 
problems  as  well  as  improving  data  locality  for  reducing  execution  time. 
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Given  below  is  a  sample  input  file  for  the  two-dimensional  problem  (3.1), 
(3.2).  Note  that  the  comments  to  the  right  of  the  numerical  values  (beginning 
with  an  “!”)  have  been  added  to  clarify  as  well  as  demonstrate  typical  input 
data.  In  actual  usage,  no  comments  will  appear  in  the  input  files.  The  run¬ 
time  options  value  2564  used  in  the  example  assures  that  all  components  of 
the  cost  will  be  calculated,  the  final  controls  with  their  spatial  locations  will 
be  saved  in  a  file,  and  the  progress  of  the  calculation  will  be  sent  to  the  output. 
The  example  in  the  next  section  (Section  3.6)  shows  how'  this  particular  value 
has  been  computed.  Also,  note  that  any  component  of  the  covariance  which 
is  entered  both  below  and  calculated  in  the  user.cov.h  include  file  will  take 
the  value  calculated  by  the  include  file,  since  such  calculations  take  place 
after  input. 


Sample  input  file  for  2D  model  (3.1)  and  cost  (3.2). 


2 

!  program  dimension 

0  16 

!  number  of  mesh  intervals 

32  32 

!  number  of  mesh  intervals 

2 

!  number  of  multigrid  sublevels 

1 

!  number  of  controls 

0 

!  non-diagonal  covariance  matrix  (0  if  false) 

0.1154700538 

!  mesh  width  h  l/\/75 

0.  0. 

!  origin  of  the  state  space  Xo  =  (0.,0.) 

2564 

!  run-time  options 

0.1154700538  0. 

!  location  of  centering  point  Xc  =  (h^O.) 

0. 

!  covariance  value  cru 

0. 

!  covariance  value  <722 

-1. 

!  control  coefficient  k(l,l)  in  drift  eqn  for  x{\ 

0. 

!  control  coefficient  k(2,l)  in  drift  eqn  for  x\2 

0.4 

!  maximum  control  value  Ui 

2 

!  number  of  cost  coefficients  c(j) 

1.  0. 

!  values  for  cost  coefficients  c{j) 

0.2  1.  1.  0.48 

\  \  H  V  a 

0.  0. 

!  underflow  cost  coefficients:  {liih] 

200.  0. 

!  overflow  cost  coefficients:  {^1,^2} 

1.  0. 

!  pi,  dimension  1  underflow  reflection 

0.  1. 

!  P2,  dimension  2  underflow  reflection 

-1.  0. 

!  9i,  dimension  1  overflow  reflection 

0.  -1. 

!  92?  dimension  2  overflow  reflection 

0. 

!  discounted  cost  factor 

150 

!  maximum  policy  updating  steps 

0.00000001 

!  stopping  criterion  tolerance 

5  5  5 

!  number  of  relaxations  per  multigrid  level 

1.2  1.2  1.2 

!  overrelaxation  parameter  for  each  level 

30 

3.6  Program  options 

At  its  simplest,  the  program  solves  the  optimal  control  problem  for  the  spec¬ 
ified  data.  This  is  option  0.  But  the  program  can  be  used  for  other  purposes. 
These  are  indicated  by  the  options  parameter  in  the  input  data.  The  user 
selects  the  options  from  the  following  list,  adds  the  values  of  the  option  indi¬ 
cator,  and  uses  this  in  the  options  input  line.  All  option  values  are  based  on 
powers  of  2  so  that  there  is  no  ambiguity  in  interpreting  the  sum  of  the  values 
of  the  individual  options.  In  all  cases  the  user  is  required  to  specify  an  option 
value  as  part  of  the  input  data.  So  to  solve  the  optimization  problem  only, 
set  the  run-time  option  to  0.  Otherwise  calculate  the  desired  option  value, 
as  described  below.  See  Table  2  for  a  summary  of  the  run-time  options. 

By  the  “components  of  the  costs,”  we  mean  the  costs  with  the  indi¬ 
vidual  Ui{-),Li{-),Ui{-),ki{-)  used  for  all  terms  with  nonzero  cost  coefficients 
Vi,  U,  Ci,  Ci,  resp.  These  costs  associated  with  the  components  are  not  weighted 
by  the  cost  coefficients,  i.e.,  the  appropriate  nonzero  cost  coeflBcients  are  set 
to  unity. 


Program  options:  Selected  at  run  time. 

•  Option  =  0  directs  the  program  to  solve  the  problem  for  the  optimum 
system  cost. 

•  Option  =  1  directs  the  program  to  save  the  final  computed  controls  as 
output  files.  These  could  subsequently  be  used  as  control  input  files 
for  another  run  where  a  cost  is  to  be  evaluated  with  a  fixed  control 
and  some  specified  cost  function.  The  number  of  files  created  equals 
the  number  of  controls.  For  each  control  Uj,  the  program  creates  an 
output  file  controlJ.data  which  contains  the  values  of  Uj  at  the  grid 
points.  See  subsection  3.7  for  a  description  of  manner  in  which  the 
grid  points  are  traversed. 

•  Option  =  2  directs  the  program  to  read  input  control  data  from  an 
existing  file  or  files.  This  option  is  used  when  we  wish  to  either  eval¬ 
uate  the  costs  for  a  given  control  or  else  start  the  policy  updating 
method  from  a  given  control.  The  input  file  or  files  must  be  named 
controlJ.data,  where  J  =  1,...,M  <  4.  The  input  control  will  be 
used  as  the  initial  control  in  the  policy  updating  scheme  unless  the  “no 
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control  updating”  option  32  is  also  specified.  This  option  is  generally 
used  for  controls  that  were  computed  and  saved  at  a  previous  run,  in 
particular  a  run  where  option  1  was  used. 

•  Option  =  4  directs  the  program  to  explicitly  compute  all  cost  compo¬ 
nents  (for  which  the  weights  are  non-zero)  for  the  problem  once  the 
optimal  control  problem  has  been  solved.  This  is  done  with  the  con¬ 
trol  fixed  at  its  value  at  the  final  iteration.  This  option  is  important 
for  getting  full  use  out  of  the  data.  Generally,  the  components  of  the 
optimal  costs  (with  weights  set  to  unity)  are  as  important  as  the  op¬ 
timal  cost  itself.  For  example,  see  [5],  where  one  would  generally  like 
to  know  what  the  actual  mean  buffer  overflow,  mean  control  loss  and 
mean  buffer  waiting  time  are  for  the  given  optimal  control  or  an  a-priori 
given  control. 

•  Option  =  8  directs  the  program  to  evaluate  only  the  cost  components 
for  the  specified  system  with  the  control  fixed.  (It  does  not  compute  an 
optimal  control.)  Combine  with  option  =  2  or  option  =  16  if  nonzero 
control  values  are  desired.  Otherwise  a  “zero”  control  will  be  used. 
The  control  used  will  not  be  updated. 

•  Option  =  16  directs  the  program  to  use  a  given  threshold  control  when 
evaluating  system  cost.  This  option  requires  additional  inputs  for  the 
threshold  levels  See  the  discussion  under  “program  inputs.”  Data  for 
this  option  is  entered  last  of  all  in  the  input  sequence.  There  is  no 
control  updating  under  this  option.  It  is  included  simplj^  because  in 
many  applications  one  wishes  to  compute  performance  under  standard 
simple  controls,  and  to  compare  the  costs  to  that  under  the  optimal. 
With  other  options,  we  allow  arbitrary  controls  to  be  used,  either  be¬ 
cause  we  wish  to  know  the  performance  under  them  or  because  we  wish 
to  start  the  policy  updates  with  some  given  “good”  control.  This  op¬ 
tion  16  facilitates  this  process  when  the  control  is  of  a  simple  threshold 
type.  Under  this  option,  the  control  is  held  fixed  through  the  program 
execution  (no  policy  updating). 

•  Option  =  32  directs  the  program  to  do  no  control  updating.  If  this 
option  is  used  alone  (or  with  option  4),  then  the  program  computes  the 
cost  for  the  zero  control.  It  also  computes  each  of  the  components  of  the 
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cost,  with  unity  weights.  When  this  option  is  used  in  conjunction  with 
option  =  2,  the  control  control  remains  fixed  as  read  from  the  input 
file  throughout  the  program’s  execution.  Note  that  when  the  control 
is  not  updated,  the  variable  “maximum  number  of  policy  updates”  is 
the  maximum  number  of  total  multigrid  cycles  allowed,  so  it  is  still  an 
important  parameter. 

•  Option  =  64  directs  the  program  to  omit  the  centering  point  in  eval¬ 
uating  the  cost  and  its  components.  This  cannot  be  elected  for  the 
ergodic  cost  function.  This  option  might  be  enabled  for  the  reflect¬ 
ing  boundary  case  with  the  discounted  cost  function.  The  absorbing 
boundary  problem  (option  128)  does  not  use  the  centering  point.  A 
sample  point  location  Xs  is  still  required  input  for  the  problem  for 
output  purposes. 

•  Option  =  128  is  used  when  the  cost  function  is  for  absorbing  boundary 

Since  there  are  no  reflection  directions  for  this  model,  the  number  of 
required  inputs  is  reduced.  No  centering  point  Xc  is  used  for  this 
model  by  default.  A  sample  point  location  Xg  is  still  required  as  input 
for  the  problem  since  the  output  is  given  as  If  option  1024  is 

also  specified  then  l/(a;)  -  V(A^s)  will  be  output.  See  Section  6  for  an 
example. 

•  Option  =  256  is  used  when  the  initial  model  is  defined  as  a  controlled 
Markov  chain,  not  as  a  reflected  diifusion  process.  See  Sections  7  and  8 
for  a  more  complete  description  of  the  specifics  for  this  option. 

•  Option  =  512  is  used  mainly  when  we  wish  to  keep  the  control  data 
for  postprocessing,  such  as  plotting.  The  option  directs  the  program  to 
create  an  output  file  of  the  final  control  data  which  includes  the  spatial 
coordinates,  as  well  as  the  associated  control  values.  The  output  file 
is  named  control,  data.  For  each  spatial  coordinate,  the  corresponding 
control  is  listed.  See  subsection  3.7  for  a  more  complete  description  of 
the  file  organization. 

•  Option  =  1024  directs  the  program  to  create  a  output  file  of  the  final 
value  (cost  or  relative  cost)  data  which  includes  the  corresponding  spa¬ 
tial  coordinates.  The  output  file  is  named  value. data.  For  each  spatial 
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coordinate,  the  corresponding  value  is  listed.  See  subsection  3.7  for  a 
more  complete  description  of  the  file  organization. 

•  Option  =  2048  informs  the  user  of  the  program’s  progress  through  the 
solution  process.  One  output  line  per  policy  step  is  produced  showing 
the  present  policy  update  cycle  number,  the  Ico  norm  and  Euclidean 
I2  norm  of  the  differences  between  the  current  cycle  and  the  last,  the 
value  of  the  system  cost  at  the  current  iteration,  and  logjg  of  this  cost. 


Option 

Function 

0 

calculate  optimal  cost  and  control 

1 

write  final  control  data 

2 

read  input  control  data 

4 

calculate  optimal  cost  and  cost  components 

8 

calculate  only  cost  components 

16 

use  threshold  control 

32 

no  control  updating 

64 

no  centering  point 

128 

absorbing  boundary  conditions  problem 

256 

user-supplied  probability  updating  code 

512 

write  final  control  and  spatial  data 

1024 

write  final  value  and  spatial  data 

2048 

report  computation  progress 

Table  2:  Program  run-time  options 

Suppose  that  a  user  wants  to  compute  the  cost  and  its  components  for 
a  given  model,  save  the  final  computed  optimal  control  for  review  and  post¬ 
processing,  and  also  survey  the  algorithm’s  progress  as  it  solves  a  model. 
Enabling  the  program’s  run-time  options  for  full  system  component  evalua¬ 
tion,  output  spatial  and  final  control  values,  and  computation  progress  will 
accomplish  these  goals.  (These  options  were  exactly  the  ones  specified  for 
use  in  the  input  data  for  example  1.)  For  these  options,  the  user  adds  the 
specified  option  values  and  assigns  the  parameter  a  value  of  2564.  That  is, 
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options  =  calculate  cost  +  compute  cost  components 
+  save  final  control  data  +  iteration  progress 
=  0  +  4  +  512  +  2048  =  2564. 

is  the  parameter  for  the  desired  options.  Thus,  the  program  allows  the  user  to 
tailor  the  execution  of  the  program  to  further  the  desired  research  goals.  The 
user,  of  course,  should  exercise  discretion  when  selecting  options.  For  exam¬ 
ple,  if  the  “no  control  updating”  and  “save  spatial  and  final  control  values” 
options  are  specified,  then  the  program  creates  the  output  file  contTol.data 
that  contains  only  zero  control  values.  In  the  situation  where  conflicting 
options  are  specified  (i.e.,  options  2  and  16),  the  program  selects  the  for¬ 
mer  option  and  ignores  the  latter.  If  input  control  files  are  specified  for  use 
(option  2)  but  do  not  exist,  the  program  will  issue  a  complaint  and  stop  exe¬ 
cution.  Appropriate  program  option  choices  allow  the  user  to  wisely  manage 
the  program  direction,  data  collection,  and  associated  system  resources. 

3.7  Program  Output 

As  the  program  is  executed,  output  is  directed  to  the  terminal.  To  save 
this  output  information,  the  user  need  only  redirect  the  output  to  a  file,  as 
illustrated  in  the  example  in  Section  3.3.  The  program  output  consists  of 
three  basic  parts;  program  type  and  options,  input  parameter  values,  and 
computed  results.  The  program  type  and  options  informs  the  user  of  the 
dimensionality  and  mesh  configuration  of  the  problem  as  well  the  selected 
options  used  by  the  program,  if  any.  All  input  parameters  are  printed  by 
the  program  to  allow  verification  of  the  input  data  by  the  user.  Finally,  the 
output  prints  the  results  of  the  calculations.  If  the  computation  progress 
option  2048  has  been  selected,  the  output  will  contain  an  active  account  of 
the  progress  of  the  multigrid  method  for  solving  the  problem.  One  key  item 
of  information  is  whether  or  not  the  program  terminated  before  achieving 
the  desired  precision.  If  the  program  seems  to  be  unstable,  decrease  the 
overrelaxation  factors.  If  it  has  not  achieved  the  desired  precision,  but  seems 
stable,  increase  the  maximum  number  of  allowed  policy  updates  and/or  the 
number  of  relaxations.  User  experience  is  the  best  guide  when  assigning  these 
parameters. 
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The  output  files  for  final  control  data,  spatial  and  control  data,  and  spatial 
and  value  data  (specified  by  program  options  2,  512,  and  1024,  resp.)  are 
written  in  a  standard  row/column  format.  For  the  output  file  options  which 
include  the  spatial  data,  each  line  of  the  file  contains  the  coordinates  of  a 
specified  point  in  the  space  and  the  value  of  the  cost  or  controls  at  that  point. 
Under  option  1024,  the  sequence  of  data  for  each  line  of  such  a  file  is  the 
spatial  position  x  ~  {xi^ . . .  ,Xd)^  and  the  the  associated  cost  value  at  that 
point.  Under  option  512,  the  value  of  the  M  controls  appears  in  lieu  of  the 
cost  value. 

The  output  is  written  beginning  with  the  least-valued  point  in  the  space 
and  the  space  is  traversed  in  the  following  standard  way.  Start  with  x  — 
(xi, . . .  ,Xd),  where  each  Xi  is  at  its  lowest  value.  Then  increase  Xi  until  the 
right  boundary  is  reached.  Then  increment  X2  by  a  unit,  return  Xi  to  its 
lowest  value  and  repeat.  Continue  until  X2  reaches  its  maximum  value.  The 
increment  x^  by  a  unit  and  repeat  and  so  on  until  all  points  are  reached. 

When  option  1  is  selected  to  save  the  control  data,  no  spatial  data  is 
written  but  the  physical  space  is  still  traversed  in  the  same  manner.  Control 
data  computed  by  other  programs  for  input  to  this  program  must  use  the 
same  data  output  scheme  if  the  assignment  of  the  input  control  is  to  be  done 
correctly. 

Furthermore,  all  input  control  files  must  have  the  total  number  of  mesh 
points  used  at  the  finest  level  for  the  problem  as  the  first  line  of  input.  Thus, 
input  control  data  files  created  by  other  programs  must  likewise  include  it. 
This  value  is  used  as  a  consistency  check  for  the  model  grid. 


4  Example:  2-D  correlated  noise  model 

We  now  give  another  illustrative  example,  where  the  Wiener  processes  are 
correlated.  The  system  equations  are 

dxi{t)  =  (-A  xi{t)  +  X2it))  dt  +  dW^t)  dt 

dx2{t)  =  (A  -  Ao)  dt  -  (^  +  +  Ao)  X2{t))  dt  -  wi(s)  dt  +  dW2{t)  dt 

(4.1) 

with  covariance  matrix 

A^i  +  ij,X2  —  (A^i  -f  IXX2) 

—  (Axi  +  HX2)  (A  —  Ao)xi  +  (/i  +  /^o  ~  ^o)x2  +  Aq 
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(4.2) 


where 


H  Ao 

Xi  —  - - - - • 

A(/xo  +  Ao)  +  Ao/i 

AAo 

X2  —  - - - 

X{fio  +  ^0)  +  Xofi 

lim^E  f  [c(l)ui(5)  +  c(2)max[0,a;2(5)  -  (4.3) 

T  1  Jo 

The  user  supplied  include  files  follow. 

•  user^drift.h: 

if(  j  .eq.  1  )  then 
drift  —  — lambda'^x(l)  +  mu*x(2) 
else 

drift  =  (lambda.-lambdaO)*x(l)  -  (mu+mu0+lambda0)*x(2) 
k  +  k(2J)^u(i,l) 

endif 

Comment:  In  this  example,  the  drift  equations  have  one  control  which 
is  denoted  by  u(i,l).  Since  there  is  no  control  variable  in  the  first 
equation  but  only  in  the  second  equation,  the  input  values  for  {k(l,l), 
k(2,l)}  would  be  {0.,  -1.}  in  this  model.  The  equations  are  dimension- 
ally  indexed  by  the  variable  j  as  required. 

•  user_cost.h: 

cost(i)  =  c(l)*u(i,l)  -f  c(2)*MAX(0.,  x(2)— B) 

Comment:  This  cost  description  utilizes  the  Fortran  “MAX”  function 
to  extract  the  nonnegative  component  of  x(2)-B.  The  reflection  cost 
terms  L  and  U  have  a  fixed  structure  so  their  weights  li,Vi  are  again 
entered  as  run  time  data. 


and 


The  cost  function  is 
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•  userJnit.h: 


xlb  =  (mu'^lambdaO)  /  (lambda*(muO  +  lambdaO)  +  lambdaO*mu) 
x2b  =  (lambda'^'lambdaO)  /  (lambda*(muO  +  lambdaO)  +  lambdaO*mu) 

Comment:  Evaluate  xlb  and  x2b  once  for  the  duration  of  the  program’s 
execution.  These  values  will  be  used  in  the  user.covar.h  file  as  part  of 
the  covariance  calculations. 

•  user_covar.h: 

cov(l,l)  =  lambda*xlb  +  mu’''x2b 
cov(l,2)  =  — (lambda*xlb  +  mu*x2b) 
cov(2,2)  =  (lambda— lambdaO)*xlb 
k  +  (mu  +  muO  —  Iambda0)*x2b  +  lambdaO 

Comment:  Compute  the  covariance  terms  from  the  specified  user  pa¬ 
rameters.  Note  that  the  term  cov(2,l)  is  not  evaluated.  The  code 
will  assign  the  value  of  cov(l,2)  to  the  term  cov(2,l)  as  it  executes  to 
assure  that  the  covariance  matrix  is  symmetrical.  The  lower  diagonal 
covariance  elements  always  receive  their  assigned  values  from  the  corre¬ 
sponding  symmetrical  upper  triangular  elements  to  assure  symmetry. 
The  character  which  appears  in  the  second  covariance  equation 
description  is  placed  in  column  6  for  Fortran  statement  continuation, 
according  to  Fortran  77  convention. 

•  user_var.h: 

real  lambda,  mu,  lambdaO,  muO,  B,  xlb,  x2b 
common  lambda,  mu,  lambdaO,  muO,  B,  xlb,  x2b 


Comment:  The  variables  which  we  have  introduced  for  this  model  are 
declared  and  typed. 

•  userJn.h: 

write(STDERR,*)  ’lambda  mu  lambdaO  muO  B  ?’ 
read(5,*)  lambda,  mu,  lambdaO,  muO,  B 
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•  user_out.h: 


write(6,*)  ’lambda  =  lambda 
write(6,*)  ’lambdaO  =  lambdaO 
write(6,*)  ’  mu  =  mu 
write(6,*)  ’  muO  =  muO 
write{6,*)  ’  B  =  B 
write(6,*)  ’  xlb  =  xlb 
write(6,*)  ’  x2b  =  x2b 

Comment:  Files  userJn.h  and  user^out.h  are  the  means  by  which  we 
enter  data  to  the  code  and  output  it  for  our  specific  model  variables. 

•  user_kbdout.h: 

write(KBD,*)  lambda,  mu,  lambdaO,  muO,  B 

Comment:  This  file  provides  the  necessary  statements  for  writing  user- 
supplied  data  to  the  default  output  file  input. prompted  by  utilizing  the 
KBD  descriptor  for  the  WRITE  statements. 

Below  an  example  input  file  for  the  two-dimensional  model  is  given.  Note 
that  the  comments  to  the  right  of  the  numerical  values  (beginning  with  an 
“!”)  have  been  added  to  clarify  as  well  as  demonstrate  typical  input  data. 
In  actual  usage,  no  comments  appear  in  input  files. 
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Sample  input  file  for  2D  correlated  noise  model 


2 

63  95 
63  127 
1 
1 
1 

0.0316227766 
0.  0. 

516 
0.  0. 

0.  0. 

0. 

0. 

-1. 

1. 

2 

1.  500. 

0.5  1.  0.1  0.2  2. 
0.  0. 

0.  0. 

1.  0. 

0.  1. 

-1.  0. 

0.  -1. 

0. 

1500 

0.00000001 
5  5 

1.1  1.1 


program  dimension 

number  of  mesh  intervals  {N^ ,  N^) 

number  of  mesh  intervals  {N2  ,  N2  ) 

number  of  multigrid  sublevels 

number  of  controls 

non-diagonal  covariance  (0  if  false) 

mesh  width  h  l/v/lOOO 

origin  of  the  state  space  Xq 

run-time  options 

location  of  centering  point  Xc 

covariance  values  for  an  and  a^ 

covariance  value  for  <722 

control  coefficient  k(l,l)  for  drift  eqn  1 

control  coefficient  k(2,l)  for  drift  eqn  2 

maximum  control  value  Uj 

number  of  cost  coefficients  c(j ) 

values  for  cost  coefficients  c(j) 

\  IX  Xo  Ho  B 

underflow  cost  coefficients:  li,l2 

overflow  cost  coefficients:  V\,V2 

Pi,  dimension  1  underflow  reflection 

P2,  dimension  2  underflow  reflection 

qi ,  dimension  1  overflow  reflection 

q2,  dimension  2  overflow  reflection 

discounted  cost  factor 

maximum  policy  updating  steps 

stopping  criterion  tolerance 

number  of  relaxations  per  multigrid  level 

overrelaxation  parameter  for  each  level 
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5  Example:  3-D  problem 

Here  is  another  example,  which  arose  in  a  study  of  a  multiplexer  with  two 
user  classes  under  heavy  traffic  [5].  The  systems  equations  are: 


dxi{t)  =  {viX2{t)  +  V2Xz{t)  -  a)dt  -  (ui(s)  +  1/2(5))  dt 

+dW2{t)  +  dW3{t)  +  dL{t)-dU{t)  (5.1) 

dxj[t)  =  — (Aj_i  +  y,j-i)xj(t)  dt  +  dWij-i{t),  j  =  2, 3. 

The  H''n(-)?  H'2(-)  and  Wzi-)  are  mutually  independent  Wiener 

processes  with  variances 

\A-[  +  fil  M  ~r  1^2  J 

ElW„.,(l)f  ■  ■'■  =  2.3. 

Aj_i  +  flj-i 

We  define  E[Wi{l)f  =  £;[W2(1)]^  +  £^[1^3(1)]^.  For  ui  and  c{i)  non-negative, 
we  use  the  stationary  cost: 


EviU(l)  +  E 


Y^^c{i)ui{s)  +  c{Z)x-i{s)ds 


The  threshold  control  option  will  be  used  for  this  problem. 


•  user_drift.h: 
if(  j  .eq.  1  )  then 

drift  =  nu(l)*x(2)  -f  nu(2)*x(3)  -  a 
k  +  k(l,l)XiT)  +  k(l,2)*u(^2) 
elseif(  j  .eq.  2  )  then 
drift  =  —  (lambda(l)  -|-  mu(l))  *  x(2) 
else 

drift  =  —  (lambda(2)  -|-  mu(2))  *  x(3) 
endif 


Comment:  Here  the  model  uses  two  controls  u(i,l)  and  u(i,2)  in  its 
drift  equations.  The  necessary  vector  index  i  appears  as  required,  as 
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Pr 


well  as  the  dimension  index  j .  The  first  drift  description  statement  is 
continued  to  the  following  line  by  the  use  of  the  character  which 
is  placed  in  column  6. 

•  user  .cost. h: 

cost(i)  =  c(l)’*'u(i,l)  +  c(2)*u(i,2)  +  c(3)*x(l) 

•  user_covar.h: 
cov(l,l)  =  (s2  +  1.0) 

(  nu(l)*b(l)*lambda(l)  /  (lambda(l)  +  mu(l)) 

+  nu(2)*b(2)*lambda(2)  /  (lambda(2)  +  mu(2))  ) 
cov(2,2)  =  2.*lambda(l)*mu(l)*b(l)  /  (lambda(l)  +  mu(l)) 
cov(3,3)  =  2.*lambda(2)*mu(2)*b(2)  /  (lambda(2)  +  mu(2)) 

Comment:  The  evaluations  for  the  covariance  terms  which  are  used 
in  this  model.  Recall  that  the  program  could  read  the  numerical  val¬ 
ues  directly  from  input  if  the  user  so  desired,  in  which  case  the  file 
user-coimr.h  would  be  left  empty. 

•  user_var.h: 

real  lambda(2),  mu(2),  nu(2),  b(2),  a,  s2 
common  lambda,  mu,  nu,  b,  a,  s2 


Comment:  We  define  and  type  the  variables  which  our  Fortran  descrip¬ 
tion  of  the  model  utilizes. 

•  userJn.h: 

write(STDERR,*)  ’lambda(l)  mu(l)  nu(l)  b(l)  ?’ 
read(5,*)  lambda(l),  mu(l),  nu(l),  b(l) 
write(STDERR,*)  ’lambda(2)  mu(2)  nu(2)  b(2)  V 
read(5,*)  lambda(2),  mu(2),  nu(2),  b(2) 
write(STDERR,*)  ’a  ?’ 
read(5,*)  a 

write(STDERR,*)  ’sigma**2  V 
read(5,*)  s2 
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•  user_out.h: 


write(6,*)  ’lambda{l,2}  =  lambda(l),  lambda(2) 
write{6,*)  ’  mu{l,2}  =  mu(l),  mu(2) 
write(6,’^)  ’  nu{l,2}  =  nu(l),  nu(2) 
write(6,*)  ’  b{l,2}  =  b(l),  b(2) 

write(6,*)  ’Heavy  traffic  constant  a  =  a 
write(6,*)  ’sigma* *2  constant  =  ’,  s2 

Comment:  Allow  the  input  of  values  for  the  introduced  variables  and 
output  their  values  for  reference  by  the  program. 

•  user_kbdout.h: 

write(KBD,*)  lambda(l),  mu(l),  nu(l),  b(l) 
write(KBD,*)  lambda(2),  mu(2),  nu(2),  b(2) 
write(KBD,*)  a 
write(KBD,*)  s2 


Comment:  This  file  provides  the  necessary  statements  for  writing  user- 
supplied  data  to  the  default  output  file  input. prompted  by  utilizing  the 
KBD  descriptor  for  the  WRITE  statements. 

Below  an  example  input  file  for  this  three-dimensional  model  is  given.  Note 
that  the  comments  to  the  right  of  the  numerical  values  (beginning  with 
an  “!”)  have  been  added  to  clarify  as  well  as  demonstrate  typical  input 
data.  In  actual  usage,  no  such  comments  appear  in  input  files.  Note  again 
that  although  the  an  are  specified,  they  may  be  overwritten  by  the  calcula¬ 
tions  given  in  the  user-cov.h  include  file.  The  reflection  directions  for  this 
model  are  normal  to  the  surface  faces. 


3 

0  16 
32  32 
32  32 
2 
2 


Sample  input  file  for  3D  multiplexer  model 

!  dimension  of  state 
!  number  of  mesh  intervals  (A^C,  A^i") 
!  number  of  mesh  intervals  {N2  ,  N2) 
!  number  of  mesh  intervals 
!  number  of  multigrid  sublevels 
!  number  of  controls 


0 

0.05 
0.  0.  0. 

564 

0.05  0.  0. 

0. 

0. 

0. 

-1.  -1. 

0.  0. 

0.  0. 

0.2  0.4 
3 

2.5  5.  1. 

0.4  2.  1.  0.5 
0.4  1.  1.  0.5 
0.48 
0. 

0.  0.  0. 

10.  0.  0. 

1.  0.  0. 

0.  1.  0. 

0.  0.  1. 

-1.  0.  0. 

0.  -1.  0. 

0.  0.  -1. 

0. 

50 

0.00000001 
3  5  5 

1.2  1.2  1.2 

1.4  2 

1.5  3 


!  ofF-diagonal  covariance  matrix  (0  if  false) 
!  mesh  width  h  =  l/y/W 
!  origin  of  the  state  space  Xq 
!  run-time  options 
location  of  centering  point  Xc 
covariance  value  a\\ 
covariance  value  (T22 
covariance  value  azz 
drift  eqn  1  control  coeffs  k(l,l),  k(l,2) 
drift  eqn  2  control  coeffs  k(2,l),  k(2,2) 
drift  eqn  3  control  coeffs  k(3,l),  k(3,2) 
maximum  control  values  u\  and  U2 
number  of  cost  coefficients  c{j) 
values  for  cost  coefficients  c{j) 

Ai  i/i  bi 

X2  Hz  ^2  62 

heavy  traffic  constant  a 
sigma’''*2  constant 
underflow  cost  coefficients:  li,l2ih 
overflow  cost  coefficients:  vi,V2,V3 
Pi  dimension  1  underflow  reflection 
P2  dimension  2  underflow  reflection 
Pz  dimension  3  underflow  reflection 
qi  dimension  1  overflow  reflection 
52  dimension  2  overflow  reflection 
qz  dimension  3  overflow  reflection 
discounted  cost  factor 
maximum  policy  updating  steps 
stopping  criterion  tolerance 
number  of  relaxations  per  multigrid  level 
overrelaxation  parameter  for  each  level 
control  1  threshold:  active  if  2:2  >  1.4 
control  2  threshold:  active  if  2:3  >  1.5 
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6  Example:  absorbing  boundary  problem 


The  following  example  of  a  two  dimensional  control  problem  with  an  ab¬ 
sorbing  boundary  will  illustrate  how  the  program  is  used  for  this  case.  The 
system  is 

dx\{t)  =  {ux2{t)  —  a)dt  —  Ui{s)dt  +  dWi  (t), 

dx2{t)  =  -(A  -h  ij)x2{t)  dt  +  dW2{t). 

The  covariance  matrix  is  diagonal  with  elements  cr(l,l)  =  0  and  <t(2,2)  = 
2  A  /i/(A  +  n).  The  (undiscounted,  3  =  0)  cost  is 


W(x,u)  =  E  f  c(l)ui(5)d5+  /  c{2)xi{s)ds  +  Eg{xr)., 
L./o  Jo 


where  r  =  min{t  :  x{t)  ^  G.  The  boundary  cost  g{x)  is  constant  on  the  set 
where  Xi  =  it  equals  x^  +  xl  where  X2  =  A"^  and  it  equals  |xi  -f  X2\ 
elsewhere. 


•  user_drift.h: 
if(  j  .eq.  1  )  then 

drift  =  nu*x(2)  —  a  +  k(l,l)'*'u(i,l) 
else 

drift  =  —  (lambda  +  mu)  *  x(2) 
endif 


Comment:  No  changes  for  the  drift  equations  from  the  earlier  multi¬ 
plexer  example. 

•  user.cost.h: 

cost(i)  =  c(l)'^u(i,l)  +  c(2)*x(l) 

Comment.  Note  that  the  part  of  the  cost  dealing  with  the  absorbing 
boundary  described  by  g{x)  is  not  included  here  because  its  structure 
is  fixed.  Otherwise,  the  cost  rate  k{x,  u{x))  is  the  same  as  we  used  in 
the  2D  multiplexer  example.  The  control  variable  includes  the  index  i 
as  required  in  all  models. 
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•  user.covar.h: 


cov(2,2)  =  2.*lambda*mu  /  (lambda  +  mu) 

Comment:  The  evaluation  of  the  nonzero  covariance  is  the  same  as  in 
the  2D  multiplexer  example. 

•  user_boundary.h 

if(  x(l)  .eq.  xmax(l)  )  then 
value(i)  =  xlbcost 
elseif(  x(2)  .eq.  xmax(2)  )  then 
value(i)  =  x(l)**2  +  x(2)**2 
else 

value(i)  =  ABS(  x(l)  +  x(2)  ) 
endif 

Comment:  We  assign  the  respective  costs  to  the  absorbing  boundary 
as  required  by  our  function  g{x).  The  spatial  values  Xj  of  each  bound¬ 
ary  point  are  available  for  use  as  well  as  the  maximum  and  minimum 
boundary  values  for  each  dimension  j  (xmax(j)  and  xmin(j),  respec¬ 
tively).  Note  the  required  use  of  the  vector  index  i  with  the  value 
array.  The  values  for  the  xmax  and  xmin  arrays  are  strictly  local  to 
the  absorbing  boundary  evaluation  routine  and  cannot  be  utilized  else¬ 
where. 

•  user.var.h: 

real  lambda,  mu,  nu,  a,  xlbcost 
common  lambda,  mu,  nu,  a,  xlbcost 

Comment:  In  addition  to  the  variables  used  for  the  drift  and  covariance 
terms  of  the  multiplexer  model,  the  variable  “xlbcost”  is  introduced 
for  assigning  a  constant  cost  to  the  Xi  overflow  boundary. 

•  userJn.h: 

write(STDERR,*)  ’lambda  mu  nu  a  ?’ 
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read(5,*)  lambda,  mu,  nu,  a 
write(STDERR,*)  ’xl  overflow  cost  constant  ?’ 
read(5,*)  xlbcost 

•  user  .out.  h: 

write(6,*)  ’lambda  =  ’,  lambda 
write(6,*)  ’mu  =  ’,  mu 
write(6,*)  ’nu  =  ’,  nu 
write(6,*)  ’a  =  ’,  a 

write(6.*)  ’xl  overflow  cost  constant  =  ’,  xlbcost 


Comment:  The  input  and  output  statements  for  this  model  again  allow 
the  input  and  output  of  user-defined  variables  needed  in  describing  the 
model. 

•  userJvbdout.h: 

write(KBD,’*')  lambda,  mu,  nu,  a 
write(KBD,*)  xlbcost 

Comment:  This  file  provides  the  necessary  statements  for  writing  user- 
supplied  data  to  the  default  output  file  input. prompted  by  utilizing  the 
KBD  descriptor  for  the  WRITE  statements. 

Given  below  is  a  sample  input  file  for  the  two-dimensional  absorbing  bound¬ 
ary  problem.  The  principle  difference  for  inputs  to  this  type  of  problem  is 
that  there  are  no  overflow  and  underflow  cost  coefficients  and  no  boundary 
reflections  terms.  In  general,  the  input  file  will  require  two  less  input  lines 
for  the  overflow  and  underflow  costs  and  (2*dim)  fewer  input  lines  for  the 
boundary  reflections.  Hence,  this  2D  absorbing  example  has  deleted  six  lines 
associated  with  the  boundaries  which  appeared  in  the  2D  reflecting  case. 
As  before,  comments  have  been  added  for  descriptive  purposes  and  do  not 
appear  in  actual  input  file  usage. 
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Sample  input  file  for  2D  absorbing  boundary  model 


2 

0  16 
32  32 
2 
1 
0 

0.1154700538 
0.  0. 

644 

0.1154700538  0. 
0. 

0. 

-1. 

0. 

0.4 

2 

1.  3. 

0.2  1.  1.  0.48 
200. 

0. 

150 

0.00000001 
5  5  5 

1.2  1.2  1.2 


program  dimension 
number  of  mesh  intervals  (A^f,  iVj'*’ ) 
number  of  mesh  intervals 
number  of  multigrid  sublevels 
number  of  controls 

non-diagonal  covariance  matrix  (0  if  false) 
mesh  width  h  ^  l/\/75 
origin  of  the  state  space  Xq  =  (0.,0.) 
run-time  options 

location  of  sample  point  Xs  —  (/i,0.) 
covariance  value  an 
covariance  value  (T22 

control  coefficient  k(l,l)  in  drift  eqn  for  a’(l) 
control  coefficient  k(2,l)  in  drift  eqn  for  x{2) 
maximum  control  value  Ui 
number  of  cost  coefficients 
values  for  cost  coefficients 
\  V  a 

absorbing  boundary  cost  for  Xi  overflow 

discounted  cost  factor 

maximum  policy  updating  steps 

stopping  criterion  tolerance 

number  of  relaxations  per  multigrid  level 

overrelaxation  parameter  for  each  level 


7  Continuous  Time  Markov  Chain  Control 
Problems 


The  previous  discussion  concerned  control  problems  for  diffusion  type  models 
such  as  (2.1).  Since  the  basic  computational  method  involves  approximation 
by  appropriate  controlled  Markov  chains,  it  is  natural  that  the  codes  can  be 
effectively  used  on  problems  that  are  posed  originally  as  control  problems  on 
Markov  chain  models.  The  basic  details  are  the  same  for  the  reflecting  and 
the  absorbing  boundary  cases.  The  main  discussion  will  be  for  the  reflect¬ 
ing  boundary  case,  and  then  the  few  required  alterations  for  the  absorbing 
boundary  will  be  given,  as  for  the  diffusion  model. 
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The  structure  of  the  state  space  is  the  same  as  discussed  for  the  Markov 
chain  approximations  in  Section  2.  It  is  an  A— grid  Gh  on  a  d—  dimensional 
hyperrectangle  G  =  with  each  point  communicating  only  to 

its  immediate  neighbors,  to  be  further  described  below.  The  value  of  h  is  an 
input  parameter.  We  assume  the  following  boundary  behavior,  analogous  to 
what  was  used  for  (2.1).  If  a  point  is  on  the  boundary  and  tries  to  move  out  of 
Gh,  then  it  is  immediately  reflected  back  with  the  mean  reflection  directions 
being  described  as  they  were  in  Section  2,  by  the  vectors  Pi,qi,i  =  1, . . .  ,d. 
Indeed,  this  behavior  fits  many  continuous  time  Markov  chain  models  which 
arise  in  telecommunications  applications. 

For  each  point  x  £  Gh,  the  behavior  is  defined  in  terms  of  controlled 
transition  rates  r{x,  y\u{x)),  where  y  varies  over  the  neighbors  of  x.  The  rates 
are  linear  in  the  controls  and  there  can  be  up  to  four  controls,  exactly  as  in 
Section  2.  Let  e,-  denote  the  unit  vector  in  the  i-th  coordinate  direction. 
Then,  for  x  £  Gh,  the  rates  for  movement  along  the  coordinate  axes  are 
wi'itten  as 

M 

r{x,x  +  eih\u)  =  rf{x)  +  '^  k^^Um,  /=l,...,d,  (7.1) 

m=l 

M 

r{x,x  —  eih\u)  rf  (x)  +  k^Um,  /  =  l,...,d,  (7-2) 

m=l 

with  analogous  formulas  for  the  “off-coordinate  axes”  rates  such  as  r{x,x  + 
o-ih  -|-  emh\u{x))  which  the  user  will  define  in  the  include  file  user jprohs.h, 
as  illustrated  below.  One  significant  difference  between  this  type  of  model 
and  the  diffusion  model  is  that  the  parameters  k^  are  not  system-defined 
variables.  Thus,  the  user  must  either  use  explicit  constants  for  all  the  cost 
coefficients  in  the  file  user  .cost. h,  or  else  define  an  array  into  which  will  be 
put  input  data  for  these  variables.  As  in  Section  2,  there  are  such  that 
0  <  ui{x)  <  ui. 

The  cost  function.  For  the  diffusion  model  case  with  a  reflecting  bound¬ 
ary,  the  cost  (2.3)  or  (2.4)  was  specified  as  the  sum  of  two  components:  a 
component  due  to  the  reflection  and  one  due  to  a  cost  rate  k(x,  u(x)).  For  the 
continuous  time  Markov  chain  problem,  the  cost  is  still  defined  as  the  sum 
of  two  components.  The  component  due  to  the  cost  rate  k(-)  is  as  in  (2.3) 
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and  (2.4).  The  cost  associated  with  boundary  reflections,  is  of  course  ac¬ 
crued  only  when  the  process  attempts  to  leave  the  state  space  and  must  be 
returned  by  a  reflection.  As  in  the  diffusion  case,  the  user  will  input  data 
for  the  overflow  and  underflow  costs,  as  well  the  reflection  directions.  The 
reflecting  boundary  cost  structure  differs  from  what  was  used  in  the  diffusion 
case  in  the  following  way:  the  inputs  li,  Ui  are  now  the  actual  cost  increments 
associated  with  a  reflection  step,  according  to  the  boundary  on  which  the 
reflection  occurs. 

Transitions  in  “off-coordinate  axes”  directions.  We  now  describe  the 
format  for  the  transition  rates  in  ‘off-coordinate  axes  directions”  such  as 
r(x,x  -h  eih  -j-  tmh\u).  I  <  m.  In  all  such  expressions  the  arguments  for  each 
vector  pair  {ei,em}  will  be  written  in  lexicographic  order,  i.e.,  I  <  m.  The 
Markov  chain  approximation  to  the  model  (2.1)  involved  such  transitions 
only  if  the  covariance  matrix  were  not  diagonal  [7].  Recall  that  one  of  the 
data  entries  told  the  program  whether  the  covariance  matrix  was  diagonal 
or  not.  Suppose  that  d  —  2,  and  ai2  >  0.  Then  the  transition  probabilities 
for  the  Markov  chain  approximation  for  the  diffusion  model  satisfy  p{x,  x  -|- 
eih  -f  e2h\u)  >  0,  p[x,x  —  e\h  —  e2h\u)  >  0,  but  for  the  transitions  for  the 
other  directions,  we  h.?ive  p{x,x-\-eih-e2h\u)  =  0,  p{x ,  x  -  tih  +  e2h\u)  =  0, 
and  conversely  if  <Ti2  <  0.  For  the  current  Markov  chain  model,  we  use  a 
similar  structure  (so  that  the  same  code  can  be  used).  Thus,  for  d  —  2, 
either  the  {r{x,x  +  eyh  fl-  e2h\u),  r{x,x  —  e-[h  —  e2h\u)]  can  be  nonzero  or 
{r(3;,  X  eih  —  e2h\u),  r{x,  x  —  tih  e2h\u)}  can  be  nonzero.  For  general  d, 
for  each  pair  /  <  m,  either  {r(a:,  x  +  eih-\-  emh\u),  r{x,  x  -  eih  -  erah\u)]  can 
be  nonzero  or  {r(a;,  x  -f  tih  —  emh\u)^  r{x,  x  —  eih  +  emh\u)'\  can  be  nonzero. 

The  choice  of  the  directions  (“northeast,  southwest”  or  “northwest,  south¬ 
east  )  is  specified  by  using  the  covariance  matrix  used  for  the  diffusion  model. 
The  covariance  itself  has  no  meaning  for  this  continuous  time  Markov  chain 
model.  But  the  signs  of  the  off-diagonal  (the  upper  triangular)  ai^rn  will 
be  used  to  indicate  the  directional  choices,  as  follows.  Thus,  if  there  are 
off-coordinate  transitions  possible  for  the  chain,  then  for  the  “nondiagonal 
covariance  matrix”  line  in  the  input  data,  write  1.  The  signs  of  the  sym¬ 
metrical  off-diagonal  aim  determine  the  allowed  directions.  The  input  values 
for  the  off-diagonal  aim-,  I  <  can  be  any  values,  provided  only  that  their 
signs  are  correct.  For  example,  if  the  directions  {e;  -f  e^,  — e;  —  e^}  are  to 
be  selected,  then  aim  >  0,  and  conversely  if  {c;  —  e^,  — e/  -f  e^}  are  to  be 
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selected  then  aim  <  0.  The  input  values  for  the  diagonal  an  are  unimportant 
so  they  can  be  any  non-negative  number. 

Defining  the  transition  rates.  The  formulas  for  the  transition  rates  are 
provided  by  the  user  in  the  file  userjprohs.h.  The  system  variables  rlft(i,j) 
and  rrht(i,j)  are  used  only  when  the  model  is  a  controlled  Markov  chain,  as 
described  in  this  section.  Evaluate  the  transition  rate  of  a  point  x  to  its 
neighbors  x  —  ejh  in  the  array  rlft(i,j)  and  to  a:;  -|-  Cjh  in  the  array  rrht(i,j)  as 
follows: 


rlft(i,j)  =  r{x,  X  —  ejh\u),  for  j  =  l, . . . ,  d 

and 

rrht(i,j)  =  r{x,  x  -f  ejh\u),  for  j=l, . . . ,  d. 

Note  that  the  index  i  is  used  by  the  system  as  before  for  indexing  the  vector 
spatial  representation  of  the  grid  data.  It  must  be  included  in  the  file  names, 
in  order  for  the  program  to  interpret  the  file  properly.  The  dimensional  index 
j  must  be  specified  by  the  user.  The  formulas  for  evaluating  the  transition 
rates  can  be  as  simple  or  complex  as  the  model  requires. 

The  program  uses  the  same  format  for  the  specification  of  the  transition 
rates  in  the  “diagonal”  directions.  Arrays  rrht(i,j)  and  rlft(i,j)  are  used  for 
these  values,  where  j=d-|-l,...,  d+D,  where  D  is  just  the  binomial  coefficient 
“choosing  2  out  of  d.”  In  detail,  the  index  j  has  the  following  meaning: 

•  Dimension  d=2: 

j  =  3  for  {61,62},  D=l, 

•  Dimension  d=3: 

j  4,5,6  for  {61,62},  {61,63},  {62,63},  respectively,  D=3, 

•  Dimension  d=4: 

j  =  5,6,7  for  {61,62},  {61,63},  {61,64},  respectively,  and 
j  =  8,9,10  for  {62,63},  {62,64},  {63,64},  respectively,  D=6. 

For  atm  >  0,  the  off-coordinate  axes  transitions  are  assigned  as 
rrht(i,j)  =  r{x,  x  -|-  eih  -|-  emh\u),  for  j=d  +  1, . . .  ,d  +  D 
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and 


rlft(i,j)  =  r{x,  X  —  eih  —  emMu),  for  j=d  +  1, . . .  ,d  +  D 
Otherwise  for  aim  ^  0,  the  off-coordinate  axes  transitions  are  assigned  as 
rrht(i,j)  =  r{x,  x  -f  eih  —  em/?|u),  for  j=d  +  1, . . .  ,d  D 

and 

rlft(i,j)  =  r{x,x  —  eih  +  emh\u),  for  j=d  +  . . .  ,d  +  D. 

8  Example:  continuous  Markov  chain 

In  this  example  we  will  use  reflecting  boundary  conditions  with  a  cost  as¬ 
signed  to  the  a’2  overflow  boundary.  The  transition  rates  are 

r{x,x  -p  e2h\u)  =  Aq  {NP  —  Xi  —  X2)  (1  —  U]) 
r{x,x  -  e2h\u)  =  ^iqX2 

/  I  I  \\  ) 

r{x,x  +  ei  —  e2n.|n)j  =  11x2 
r{x,x  —  Cl  -f  e2h\u)  =  A 

and  the  cost  function  is  similiar  to  the  correlated  noise  example  but  with  the 
addition  of  the  X2  boundary  overflow  cost  term: 

Et')2Uil) f  [c(l)ui(s) -|- c(2) max[0, 0:2(5)  —  5]] ds. 

T  T  Jo 

Recall  that  the  state  space  of  the  Markov  chain  is  an  /i— grid. 

•  user-drift. h: 

•  user_covar.h: 

Comment:  These  files  are  left  empty  since  they  are  unused  by  the  code 
for  the  continuous  time  Markov  chain  model.  The  files  must  be  present 
for  the  compilation  process  to  succeed. 
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•  user.probs.h: 

rrht(i,2)  =  lambdaO  *  (NP-x(l)-x(2))  *  (  1.0  -  u(i,l)  ) 
rlft(i,2)  =  muO  *  x(2) 
rrht(i,3)  =  mu  *  x(2) 
rlft(i,3)  =  lambda  *  x(l) 

Comment:  The  transition  rates  for  rrht(i,2)  and  rlft(i,2)  are  easily 
coded.  Since  there  are  no  transitions  along  the  first  coordinate  direc¬ 
tion,  we  omit  these  zero-value  assignments  for  efficiency.  A  covariance 
value  o'i2  ^  0.  is  needed  since  the  ofF-coordinate  axes  transitions  are 
in  the  southeast  and  northwest  directions  for  rrht(i,3)  and  rlft(i,3),  le- 
spectively. 

•  user -Cost,  h: 

cost(i)  =  c(l)*u(i,l)  +  c(2)*max(0.,  x(2)-B) 

•  user.var.h: 

real  lambda,  mu,  lambdaO,  muO,  B 
integer  NP 

common  lambda,  mu,  lambdaO,  muO,  B,  NP 

•  userJn.h: 

write(STDERR,*)  ’lambda  mu  lambdaO  muO  B  NP  ?’ 
read(5,*)  lambda,  mu,  lambdaO,  muO,  B,  NP 

•  user  .out.  h: 

write(6,’^)  ’lambda  =  ’,  lambda 
write(6,*)  ’lambdaO  =  ’,  lambdaO 
write(6,*)  ’  mu  =  ’,  mu 
write(6,*)  ’  muO  =  ’,  muO 
write(6,*)  ’  B  =  ’,  B 
write(6,*)  ’  NP  =  ’,  NP 
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Comment:  As  in  previous  examples,  we  define  the  model  variables  and 
the  necessary  input  and  output  statements  so  that  our  specific  data, 
will  be  read  and  subsequentially  written  for  reference. 

•  user_kbdout.h: 

write(KBD,*)  lambda,  mu,  lambdaO,  muO,  B,  NP 


Comment:  This  file  provides  the  necessary  statements  for  writing  user- 
supplied  data  to  the  default  output  file  input. prompted  by  utilizing  the 
KBD  descriptor  for  the  WRITE  statements. 

Again,  it  is  important  to  note  that  the  control  coefficient  matrix  K  with 
elements  k(i,j)  is  not  used  for  the  continuous  time  Markov  chain  model  as  it 
was  in  the  the  diffusion  models.  Thus,  input  files  for  these  models  will  not 
have  data  for  the  k(i,j)  control  coefficients. 

Sample  input  file  for  2D  continuous  time  Markov  chain  model 


2 

!  program  dimension 

0  50 

!  number  of  mesh  intervals  (Aq“,  ) 

0  60 

!  number  of  mesh  intervals  (W^, 

1 

!  number  of  multigrid  sublevels 

1 

!  number  of  controls 

1 

!  off-diagonal  covariance  terms  (0  if  false) 

0.1 

!  mesh  width  h  —  l/-\/l00 

0.  2. 

!  origin  of  the  state  space  A^o 

772 

!  run-time  options 

4.  7. 

!  location  of  centering  point  Xc 

0.  0. 

!  covariance  values  for  a\i  and  (7\2 

0. 

!  covariance  value  for  022 

0.5 

!  maximum  control  value  Ui 

2 

!  number  of  cost  coefficients  c{j) 

15.  2. 

!  values  for  cost  coefficients  c{j) 

.5  1.  .1  .2  6.3  200 

\  \  p.  Xq  p.Q  B  NP 

0.  0. 

!  underflow  cost  coefficients:  li,l2 

0.  1. 

!  overflow  cost  coefficients:  Vi,V2 

1.  0. 

!  pi,  dimension  1  underflow  reflection 

0.  1. 

!  p2,  dimension  2  underflow  reflection 

-1.  0. 

!  9i,  dimension  1  overflow  reflection 

0.  -1. 

!  92,  dimension  2  overflow  reflection 

54 
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.00000001 
5  5 

1.0  1.0 


!  discounted  cost  factor 
!  maximum  policy  updating  steps 
!  stopping  criterion  tolerance 
!  number  of  relaxations  per  multigrid  level 
!  overrelaxation  parameter  for  each  level 


Absorbing  boundary.  The  absorbing  boundary  for  a  continuous  Markov 
chain  follows  exactly  the  same  format  as  the  diffusion  model  case.  One  need 
not  specify  the  reflection  directions  p*  and  qi  and  reflection  costs  /,■  and  Vi 
but  the  file  for  the  boundary  cost  g{x)  must  be  provided  as  before. 
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